% 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
Best Answer