%-------------------------
% 定步長的龍格庫塔4階算法
%-------------------------
h=0.02;
y1(1)=1,y2(1)=0,y3(1)=-1;
m=0;
for t=[h:h:2]
    m=m+1;
    k11=-21*y1(m)+19*y2(m)-20*y3(m);
    k21= 19*y1(m)-21*y2(m)+20*y3(m);
    k31= 40*y1(m)-40*y2(m)-40*y3(m);
    
    k12=-21*(y1(m)+h/2*k11)+19*(y2(m)+h/2*k21)-20*(y3(m)+h/2*k31); 
    k22=19*(y1(m)+h/2*k11)-21*(y2(m)+h/2*k21)+20*(y3(m)+h/2*k31);
    k32=40*(y1(m)+h/2*k11)-40*(y2(m)+h/2*k21)-40*(y3(m)+h/2*k31);
    
    k13=-21*(y1(m)+h/2*k12)+19*(y2(m)+h/2*k22)-20*(y3(m)+h/2*k32);
    k23= 19*(y1(m)+h/2*k12)-21*(y2(m)+h/2*k22)+20*(y3(m)+h/2*k32);
    k33=40*(y1(m)+h/2*k12)-40*(y2(m)+h/2*k22)-40*(y3(m)+h/2*k32);
    
    k14=-21*(y1(m)+h*k13)+19*(y2(m)+h*k23)-20*(y3(m)+h*k33);
    k24=19*(y1(m)+h*k13)-21*(y2(m)+h*k23)+20*(y3(m)+h*k33);
    k34=40*(y1(m)+h*k13)-40*(y2(m)+h*k23)-40*(y3(m)+h*k33);
    
    y1(m+1)= y1(m)+(h/6)*( k11 + 2*k12 + 2*k13 +k14);
       
    y2(m+1)= y2(m)+(h/6)*( k21 + 2*k22 + 2*k23 +k24);
     
    y3(m+1)= y3(m)+(h/6)*( k31 + 2*k32 + 2*k33 +k34);
   
end 
t=[0 : h : 2];

%----------
% 解析解
%----------
x1=1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)+1/2*exp(-40*t).*sin(40*t) ;
x2=-1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)-1/2*exp(-40*t).*sin(40*t);
x3=exp(-40*t).*sin(40*t)-exp(-40*t).*cos(40*t);


%--------------------------
% 解析解與數值解差值比較 
%--------------------------

subplot(3,1,1)
plot(t,y1-x1,'-*'),xlabel('t'),ylabel('y1')          

subplot(3,1,2)
plot(t,y2-x2,'-*'),xlabel('t'),ylabel('y2')

subplot(3,1,3)
plot(t,y3-x3,'-*'),xlabel('t'),ylabel('y3')