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
