%
%  Math 551 5/9/13
%
% demos with ode23tx, etc
% from sec. 7.7 and 7.9
%
%
F=@(t,y) 0

F = 

    @(t,y)0

ode23tx(F,[0 10],1)
[tout,yout]=ode23tx(F,[0 10],1)

tout =

     0
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10


yout =

     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1

F=@(t,y) y

F = 

    @(t,y)y

[tout,yout]=ode23tx(F,[0 2],1)

tout =

         0
    0.0800
    0.2800
    0.4800
    0.6800
    0.8800
    1.0800
    1.2800
    1.4800
    1.6800
    1.8800
    2.0000


yout =

    1.0000
    1.0833
    1.3231
    1.6159
    1.9735
    2.4103
    2.9438
    3.5954
    4.3912
    5.3631
    6.5501
    7.3852

ode23tx(F,[0 2],1)
exp(2)

ans =

    7.3891

F=@(t,y) y^2

F = 

    @(t,y)y^2

ode23tx(F,[0 2],1)
{Warning: Step size 3.512299e-15 too small at t = 1.001616e+00.} 

% 2 x 2 system: harmonic oscillator
F=@(t,y) [y(2); -y(1)]

F = 

    @(t,y)[y(2);-y(1)]

ode23tx(F,[0 2*pi],[1;0])
[tout,yout]=ode23tx(F,[0 2*pi],[1;0])

tout =

         0
    0.0001
    0.0005
    0.0025
    0.0125
    0.0625
    0.1778
    0.3416
    0.5466
    0.7911
    1.0797
    1.3461
    1.5208
    1.6963
    1.8414
    2.0304
    2.2589
    2.5288
    2.8111
    3.0544
    3.1935
    3.3219
    3.4864
    3.6922
    3.9375
    4.2270
    4.4926
    4.6660
    4.8401
    4.9862
    5.1761
    5.4055
    5.6762
    5.9577
    6.2002
    6.2832


yout =

    1.0000         0
    1.0000   -0.0001
    1.0000   -0.0005
    1.0000   -0.0025
    0.9999   -0.0125
    0.9980   -0.0624
    0.9842   -0.1769
    0.9422   -0.3350
    0.8542   -0.5197
    0.7028   -0.7110
    0.4713   -0.8814
    0.2225   -0.9742
    0.0497   -0.9980
   -0.1252   -0.9913
   -0.2672   -0.9628
   -0.4433   -0.8954
   -0.6346   -0.7715
   -0.8172   -0.5743
   -0.9446   -0.3238
   -0.9946   -0.0866
   -0.9970    0.0522
   -0.9821    0.1794
   -0.9395    0.3378
   -0.8506    0.5226
   -0.6981    0.7134
   -0.4651    0.8828
   -0.2170    0.9737
   -0.0457    0.9965
    0.1275    0.9894
    0.2702    0.9602
    0.4466    0.8919
    0.6376    0.7669
    0.8194    0.5683
    0.9447    0.3182
    0.9934    0.0820
    0.9967   -0.0007

plot(yout(:,1),yout(:,2),'-o')
axis equal
axis([-1.2 1.2 -1.2 1.2])
plot(yout(:,1),yout(:,2),'-o')
axis([-1.2 1.2 -1.2 1.2])
axis equal

% add some resistence
F=@(t,y) [y(2); -y(1)-0.1*y(2)]

F = 

    @(t,y)[y(2);-y(1)-0.1*y(2)]

ode23tx(F,[0 2*pi],[1;0])
[tout,yout]=ode23tx(F,[0 2*pi],[1;0]);
plot(yout(:,1),yout(:,2),'-o')
axis([-1.2 1.2 -1.2 1.2])
axis equal
clear
% sec. 7.9 stiffness
help odephas



% sec. 7.9 stiffness
delta=0.01;
F=@(t,y) y^2 - y^3;

opts = odeset('RelTol',1.e-4);
ode45(F,[0 2/delta],delta,opts);
delta=0.00001;
ode45(F,[0 2/delta],delta,opts);
ode23s(F,[0 2/delta],delta,opts);
diary off