%四階龍格庫塔方法
function[T,U]=kutta(a,b,h,c)
%計算點的個數
n=round((b-a)/h)+1;
%定義T,U的維數
t=zeros(n,1);
u=zeros(n,1);
%微分方程的初始條件
t(1)=a;
u(1)=c;
%顯示微分方程的初始條件
str=sprintf('t0=%g, u0=%g\n',t(1),u(1));
disp(str);
%對每個點用四階龍格庫塔方法計算
for i=1:n-1
    %系數
    k1=f(t(i),u(i));
    k2=f(t(i)+0.5*h,u(i)+0.5*h*k1);
    k3=f(t(i)+0.5*h,u(i)+0.5*h*k2);
    k4=f(t(i)+h,u(i)+h*k3);
    %龍格庫塔公式
    u(i+1)=u(i)+h*(k1+2*k2+2*k3+k4)/6;
    t(i+1)=t(i)+h;
    %顯示計算結果與誤差
    str=sprintf('%d t=%g u=%g orig=%g e=%g',i,t(i+1),u(i+1),orig(t(i+1)),abs(u(i+1)-orig(t(i+1))));
     if rem(i,8)==0;
         disp(str);
     end
 end
 %返回計算結果
 T=t;
 U=u;