Error using rk4sys (line 2) at least 4 inout argument required

% My Code

function [tp,yp] = rk4sys(dydt,tspan,y0,h,varargin) if nargin<4, error('at least 4 inout argument required'), end if any(diff(tspan)<=O),error('tspan not ascending order'), end n = length (tspan); ti = tspan(l);tf = tspan(n); if n == 2 t = (ti:h:tf)'; n = length(t); if t(n)<tf t(n+l) = tf;n=n+l; end else t = tspan; end tt = ti; y(l,:) = y0; np=1; tp(np) = tt; yp(np,:)= y(1,:); i=l; while(1) tend = t(np+l);hh = t(np+l) - t(np); if hh>h,hh = h; end while(l) if tt+hh>tend,hh = tend-tt; end kl = dydt (tt, y(i,:) ,varargin{:})'; ymid = y(i,:) + kl.*hh./2; k2 = dydt(tt+hh/2,ymid,varargin{:})'; ymid = y(i,:) + k2*hh/2; k3 = dydt(tt+hh/2,ymid,varargin{:})'; yend = y(i,:) + k3*hhi k4 = dydt(tt+hh,yend,varargin{:})'; phi = (kl+2*(k2+k3)+k4)/6i y(i+l,:) = y(i,:) + phi*hh; tt = tt+hhi i=i+l; if tt>=tend,break,end end np = np+l; tp(np) = tt; yp(np,:) = y(i,:); if tt>=tf,break,end end pa = 2560; h = 5; tspan = [1950 2000];

[t,p1] = rk4sys(@fun,tspan,pa,h); p = [2560,2780,3040,3350,3710,4090,4450,4850,5280,5690,6080]; plot(t,p,'*',t,p1); xlabel('t');ylabel('p'); legend('Actual','RK4') end

my function function f = fun(t,p) f= 0.0077*p; end

Answer 1

[Jan]

How do you call this code?

You have posted it in one block. (By the way: Use the "{} Code" button to apply a proper formatting - posting readable code is a big advantage in the forum...) Maybe you have written all the code in one file and start it by the "Run" button in the editor?

Better:

function main
pa = 2560;
h = 5;
tspan = [1950 2000];

[t,p1] = rk4sys(@fun,tspan,pa,h);
p = [2560,2780,3040,3350,3710,4090,4450,4850,5280,5690,6080];
plot(t,p,'*',t,p1);
xlabel('t');ylabel('p');
legend('Actual','RK4')
end

function [tp,yp] = rk4sys(dydt,tspan,y0,h,varargin)
if nargin<4, error('at least 4 inout argument required'), end
if any(diff(tspan)<=O),error('tspan not ascending order'), end
n = length (tspan);
ti = tspan(l);tf = tspan(n);
if n == 2
t = (ti:h:tf)'; n = length(t);
if t(n)<tf
t(n+l) = tf;n=n+l;
end
else
t = tspan;
end
tt = ti; y(l,:) = y0;
np=1; tp(np) = tt; yp(np,:)= y(1,:);
i=l;
while(1)  
tend = t(np+l);hh = t(np+l) - t(np);
if hh>h,hh = h; end 
while(l)
if tt+hh>tend,hh = tend-tt; end
kl = dydt (tt, y(i,:) ,varargin{:})';
ymid = y(i,:) + kl.*hh./2;
k2 = dydt(tt+hh/2,ymid,varargin{:})';
ymid = y(i,:) + k2*hh/2;
k3 = dydt(tt+hh/2,ymid,varargin{:})';
yend = y(i,:) + k3*hhi
k4 = dydt(tt+hh,yend,varargin{:})';
phi = (kl+2*(k2+k3)+k4)/6i
y(i+l,:) = y(i,:) + phi*hh;
tt = tt+hhi
i=i+l;
if tt>=tend,break,end
end
np = np+l; tp(np) = tt; yp(np,:) = y(i,:);
if tt>=tf,break,end
end

my function 
function f = fun(t,p)
f = 0.0077*p;
end