Tuesday, April 20, 2010

Runge Kutta 4th order code

function u = RK4(h,tf,x1,y1)
%to solve the ode for radiation using rk fourth order
I=tf/h;
fh = @ (x , y) (-2.2067)*(10^-12)*(y^4-(81*10^8));
u = y1;
for i = 1:I
k1=h*fh(x1+(i-1)*h,u);

k2=h*fh(x1+(i-1)*h+(h/2),u+(k1/2));

k3=h*fh(x1+(i-1)*h+(h/2),u+(k2/2));

k4=h*fh(x1+(i-1)*h + h,u + k3);

u= u + (1/6*(k1 + 2*(k2 + k3) + k4));

end
u
end

No comments: