format long
%% use Runge-Kutta to solve IVP
%% y'=1+y^2
%% y(0)=0
%% 
%% solution is y=tan(x)
%% this blows up at x=pi/2=1.570796326794897...
%%

%% first, see what happens if we try to find
%% y(1.6)
[x y]=RungeKuttaMethod(0,1.6,10,0);
plot(x,y)
y(11)
[x y]=RungeKuttaMethod(0,1.6,20,0);
plot(x,y)
y(21)
[x y]=RungeKuttaMethod(0,1.6,40,0);
plot(x,y)
y(41)
[x y]=RungeKuttaMethod(0,1.6,1000,0);
plot(x,y)
y(1001)
%% we see that the solution seems to be blowing up before x=1.6
%%
%% lets see what happens if we try to find y(1.55)
[x y]=RungeKuttaMethod(0,1.55,100,0);
plot(x,y)
[x y]=RungeKuttaMethod(0,1.55,1000,0);
plot(x,y)
[x y]=RungeKuttaMethod(0,1.55,10000,0);
plot(x,y)
y(10001)
[x y]=RungeKuttaMethod(0,1.55,20000,0);
y(20001)
[x y]=RungeKuttaMethod(0,1.55,40000,0);
y(40001)
%% so it looks like the solution is stable up to x=1.55
%%
%% lets try to find an accurate approximation for where the blow up occurs
%% first, get an accurate y value for a stable point near the blow up, say x=1.55, and the use these as new IC
[x y]=RungeKuttaMethod(0,1.55,100000,0);
x0=1.55;
y0=y(100001);
%% then rerun Runge-Kutta with these new IC
%% do it several times with different widths for the x-interval
%% try to narrow down where the blow up occurs
%%
%% for example, as a sanity check, we know taking the x-interval to have a width of .05 should result in a blow up
width=.05;
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
%% try smaller widths
width=.01;
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% seems stable so increase width
width=.02;
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% still seems stable so increase width
width=.03;
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
%% solution has blown up
%% so we know the singularity (the x-value where the blow up occurs) is between 1.57 and 1.58
%%
%% lets repeat the above process but with IC starting at x0=1.57
%% first need to get an accurate y0
[x y]=RungeKuttaMethod(x0,x0+width,100000,y0);
x0=x0+width
y0=y(100001)
%% now repeat the above procedure
width=.001
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101);
%% blows up, make width smaller
width=.0001
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0002
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0003
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0004
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0005
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0006
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0007
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
[x y]=RungeKuttaMethod(x0,x0+width,1000,y0);
y(1001)
%% stable, increase width
width=.0008
[x y]=RungeKuttaMethod(x0,x0+width,100,y0);
y(101)
%% blows up!
%% so we know the singularity is between
%% 1.5707 and 1.5708
%% lets double check out work
%% i.e. generate some evidence that this is correct
[x y]=RungeKuttaMethod(0,1.55,1000,0);
y(1001)
[x y]=RungeKuttaMethod(0,1.55,10000,0);
y(10001)
[x y]=RungeKuttaMethod(0,1.55,50000,0);
y(50001)
[x y]=RungeKuttaMethod(0,1.55,100000,0);
y(100001)
%% stable, update IC
x0=1.55
y0=y(100001)
[x y]=RungeKuttaMethod(x0,1.5507,1000,y0);
y(1001)
[x y]=RungeKuttaMethod(x0,1.5507,10000,y0);
y(10001)
[x y]=RungeKuttaMethod(x0,1.5507,100000,y0);
y(100001)
%% so stable at 1.5507
[x y]=RungeKuttaMethod(x0,1.5508,1000,y0);
y(1001)
[x y]=RungeKuttaMethod(x0,1.5508,10000,y0);
y(10001)
[x y]=RungeKuttaMethod(x0,1.5508,100000,y0);
y(100001)
%% blows up at 1.5508

%%
%% why did we go through the trouble of updating IC?
%% suppose we had wanted to get things more accurate
%% say we wanted to show that the singularity is between
%% 1.5707963 and 1.5707964
%% if we started at x0=0
%% it would take too many steps to get accurate approximations at 1.5707963 and 1.5707964
[x y]=RungeKuttaMethod(0,1.5707963,1000,0);
y(1001);
[x y]=RungeKuttaMethod(0,1.5707963,10000,0);
y(10001);
[x y]=RungeKuttaMethod(0,1.5707963,100000,0);
y(100001);
%% not clear if this is stabilizing and taking steps up to 1000000 would take long to run
[x y]=RungeKuttaMethod(0,1.5707964,1000,0);
y(1001);
[x y]=RungeKuttaMethod(0,1.5707964,10000,0);
y(10001);
[x y]=RungeKuttaMethod(0,1.5707964,100000,0);
y(100001);
%% again not clear if this is blowing up or not and would take a long time to run with more steps
%% the point is -- the extra steps are only needed near the singularity!! away from the singularity, you can get more accuracy with fewer steps







