% Harj. 4 teht. 4 Matlab-istunto
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>> close all; clear
global m1 m2
m1=1; m2=1;
R01=[1;0]; R02=[0;1];
V01=[1;0]; V02=-V01;

Y0=[R01;V01;R02;V02];

%addpath c:\usr\heikki\numsym04\H
addpath /home/apiola/opetus/numsym/04/H  

kakskpl(27,Y0)

[T,Y]=ode45(@kakskpl,[0 100],Y0);
R1=Y(:,1:2);
R2=Y(:,5:6);
plot(R1(:,1),R1(:,2))
hold on
plot(R2(:,1),R2(:,2),'r'); shg;
title(['Alkunopeudet [',num2str(V01(1)),'  ',num2str(V01(2)),'] ja ['...
  ,num2str(V02(1)),'  ',num2str(V02(2)),']'])

% Uudet alkunopeudet
figure(2)
V01=[-1; 0]; V02=[1;0];Y0=[R01;V01;R02;V02];
clf
[T,Y]=ode45(@kakskpl,[0 100],Y0);
R1=Y(:,1:2);R2=Y(:,5:6); plot(R1(:,1),R1(:,2)); hold on;
plot(R2(:,1),R2(:,2),'r'); shg
title(['Alkunopeudet [',num2str(V01(1)),'  ',num2str(V01(2)),'] ja ['...
  ,num2str(V02(1)),'  ',num2str(V02(2)),']'])


figure(3)
c=0.5; V01=c*[-1; 0]; V02=c*[1;0];Y0=[R01;V01;R02;V02];
clf
[T,Y]=ode45(@kakskpl,[0 100],Y0);
R1=Y(:,1:2);R2=Y(:,5:6); plot(R1(:,1),R1(:,2)); hold on;
plot(R2(:,1),R2(:,2),'r'); shg
title(['Alkunopeudet [',num2str(V01(1)),'  ',num2str(V01(2)),'] ja ['...
  ,num2str(V02(1)),'  ',num2str(V02(2)),']'])

figure(4)
c1=0.5; V01=c1*[-1; 0]; c2=0.25; V02=c2*[1;0];Y0=[R01;V01;R02;V02];
[T,Y]=ode45(@kakskpl,[0 100],Y0);
R1=Y(:,1:2);R2=Y(:,5:6); plot(R1(:,1),R1(:,2)); hold on;
plot(R2(:,1),R2(:,2),'r'); shg
title(['Alkunopeudet [',num2str(V01(1)),'  ',num2str(V01(2)),'] ja ['...
  ,num2str(V02(1)),'  ',num2str(V02(2)),']'])


ans =

   1.00000000000000
                  0
  -0.35355339059327
   0.35355339059327
  -1.00000000000000
                  0
   0.35355339059327
  -0.35355339059327

>> 
>> 
>> 
>> 
>> f=inline('1/y^2-y*t','t','y')
[T,Y]=ode45(f,[1 2],1);
Y(end,:)

f =

     Inline function:
     f(t,y) = 1/y^2-y*t


ans =

   0.82356797018734

>> opt=odeset('reltol',1e-10,'abstol',1e-10)
[T,Y]=ode45(f,[1 2],1,opt);
Y(end,:)

opt = 

              AbsTol: 1.000000000000000e-10
                 BDF: []
              Events: []
         InitialStep: []
            Jacobian: []
           JConstant: []
            JPattern: []
                Mass: []
        MassConstant: []
        MassSingular: []
            MaxOrder: []
             MaxStep: []
         NormControl: []
           OutputFcn: []
           OutputSel: []
              Refine: []
              RelTol: 1.000000000000000e-10
               Stats: []
          Vectorized: []
    MStateDependence: []
           MvPattern: []
        InitialSlope: []


ans =

   0.82356789707243

>> [T,Y]=ode113(f,[1 2],1,opt);
>> Y(end,:)

ans =

   0.82356789696604