restart;

Here we set the number of digits that we will require correct, but this comes at a cost as it requires maple to use

more registers for storage of each result

Digits:=15;

NiM+SSdEaWdpdHNHNiIiIzo=

First we will set the order of the method

methorder:=6;

NiM+SSptZXRob3JkZXJHNiIiIic=

n:=methorder-1;

NiM+SSJuRzYiIiIm

Next we set the end points for the interval we wish to approximate the solution on where

a is the starting value where the initial condition will be set and b is the ending value.

a:=0;

NiM+SSJhRzYiIiIh

b:=1;

NiM+SSJiRzYiIiIi

The variable loc will be a running variable that will point to the value of the independent variable

t as we march across the interval from a to b.

loc:=a;

NiM+SSRsb2NHNiIiIiE=

The variable xstart is the starting for x given in our initial condidtion.

xstart:=-1;

NiM+SSd4c3RhcnRHNiIhIiI=

For teaching purposes since I am not concerned with storage I will store the approximation to the

solution and and each of the derivaties at each step to do this it seems reasonable to use a two

dimensional array of the form x[i,j] where i stands for the i^th derivative and j stands for the

step number. So x[0,0] is the value for x, or the 0^th derivative of x at step 0, or the initial condition.

x[0,0]:=xstart;

NiM+JkkieEc2IjYkIiIhRighIiI=

Here we set the number of steps

steps:=10;

NiM+SSZzdGVwc0c2IiIjNQ==

The variable h is the step size

h:=evalf((b-a)/steps);

NiM+SSJoRzYiJCIwKysrKysrKyIhIzo=

Now we set the right hand side. Since I want to let maple take the derivatives I will substitute w for x

in my equations.

****************************

g:=1+w(t)^2-t^3;

NiM+SSJnRzYiLCgiIiJGJyokLUkid0dGJTYjSSJ0R0YlIiIjRicqJEYsIiIkISIi

****************************

This sets the 0th derivative of of the right hand side to g

df[0]:=g;

NiM+JkkjZGZHNiI2IyIiISwoIiIiRioqJC1JIndHRiY2I0kidEdGJiIiI0YqKiRGLyIiJCEiIg==

This loop takes the required derivatives of the right hand side

for i from 1 to n by 1 do

df[i]:=diff(g,t$i);

od;

NiM+JkkjZGZHNiI2IyIiIiwmKiYtSSJ3R0YmNiNJInRHRiZGKC1JJWRpZmZHSSpwcm90ZWN0ZWRHRjE2JEYrRi5GKCIiIyokRi5GMyEiJA==

NiM+JkkjZGZHNiI2IyIiIywoKiQtSSVkaWZmR0kqcHJvdGVjdGVkR0YtNiQtSSJ3R0YmNiNJInRHRiZGMkYoRigqJkYvIiIiLUYsNiRGLy1JIiRHRi02JEYyRihGNEYoRjIhIic=

NiM+JkkjZGZHNiI2IyIiJCwoKiYtSSVkaWZmR0kqcHJvdGVjdGVkR0YtNiQtSSJ3R0YmNiNJInRHRiZGMiIiIi1GLDYkRi8tSSIkR0YtNiRGMiIiI0YzIiInKiZGL0YzLUYsNiRGLy1GNzYkRjJGKEYzRjkhIidGMw==

NiM+JkkjZGZHNiI2IyIiJSwoKiQtSSVkaWZmR0kqcHJvdGVjdGVkR0YtNiQtSSJ3R0YmNiNJInRHRiYtSSIkR0YtNiRGMiIiI0Y2IiInKiYtRiw2JEYvRjIiIiItRiw2JEYvLUY0NiRGMiIiJEY7IiIpKiZGL0Y7LUYsNiRGLy1GNDYkRjJGKEY7RjY=

NiM+JkkjZGZHNiI2IyIiJiwoKiYtSSVkaWZmR0kqcHJvdGVjdGVkR0YtNiQtSSJ3R0YmNiNJInRHRiYtSSIkR0YtNiRGMiIiIyIiIi1GLDYkRi8tRjQ2JEYyIiIkRjciIz8qJi1GLDYkRi9GMkY3LUYsNiRGLy1GNDYkRjIiIiVGNyIjNSomRi9GNy1GLDYkRi8tRjQ2JEYyRihGN0Y2

Here we create essentially functions of the derivatives of the right hand side

for i from 0 by 1 to n do

f[i]:=unapply((subs([w(t)=y,seq(diff(w(t),t$k)=y[k],k=1..i)],df[i])),[t,y,seq(y[k],k=1..i)]);

od;

NiM+JkkiZkc2IjYjIiIhZio2JEkidEdGJkkieUdGJkYmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCgiIiJGMSokOSUiIiNGMSokOSQiIiQhIiJGJkYmRiY=

NiM+JkkiZkc2IjYjIiIiZio2JUkidEdGJkkieUdGJkkkeV8xR0YmRiY2JEkpb3BlcmF0b3JHRiZJJmFycm93R0YmRiYsJiomOSVGKDkmRigiIiMqJDkkRjUhIiRGJkYmRiY=

NiM+JkkiZkc2IjYjIiIjZio2JkkidEdGJkkieUdGJkkkeV8xR0YmSSR5XzJHRiZGJjYkSSlvcGVyYXRvckdGJkkmYXJyb3dHRiZGJiwoKiQ5JkYoRigqJjklIiIiOSdGN0YoOSQhIidGJkYmRiY=

NiM+JkkiZkc2IjYjIiIkZio2J0kidEdGJkkieUdGJkkkeV8xR0YmSSR5XzJHRiZJJHlfM0dGJkYmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCgqJjkmIiIiOSdGNiIiJyomOSVGNjkoRjYiIiMhIidGNkYmRiZGJg==

NiM+JkkiZkc2IjYjIiIlZio2KEkidEdGJkkieUdGJkkkeV8xR0YmSSR5XzJHRiZJJHlfM0dGJkkkeV80R0YmRiY2JEkpb3BlcmF0b3JHRiZJJmFycm93R0YmRiYsKCokOSciIiMiIicqJjkmIiIiOShGOyIiKSomOSVGOzkpRjtGN0YmRiZGJg==

NiM+JkkiZkc2IjYjIiImZio2KUkidEdGJkkieUdGJkkkeV8xR0YmSSR5XzJHRiZJJHlfM0dGJkkkeV80R0YmSSR5XzVHRiZGJjYkSSlvcGVyYXRvckdGJkkmYXJyb3dHRiZGJiwoKiY5JyIiIjkoRjgiIz8qJjkmRjg5KUY4IiM1KiY5JUY4OSpGOCIiI0YmRiZGJg==

This double loop is where we actually do the TS method and march across the interval

for j from 0 by 1 to steps-1 do

for i from 1 by 1 to n do

x[i,j]:=f[i-1](loc,seq(x[k,j],k=0..i-1));

od;

x[0,j+1]:=sum(x[k,j]*h^k/(k!),k=0..n);

loc:=loc+h;

od:

Now lets have maple solve the inital value problem

deq:=diff(w(t),t)=g;

NiM+SSRkZXFHNiIvLUklZGlmZkdJKnByb3RlY3RlZEdGKTYkLUkid0dGJTYjSSJ0R0YlRi4sKCIiIkYwKiRGKyIiI0YwKiRGLiIiJCEiIg==

ic:=w(a)=xstart;

NiM+SSNpY0c2Ii8tSSJ3R0YlNiMiIiEhIiI=

Let's have the power of a thousand mathematicians try to find an analytic solution if nothing comes up

it means they gave up

dsolve({deq,ic},w(t));

So now lets have them find their numerical solution, it chooses the step size to gain an accurate solution, and in general it is pretty good

for the problems here we will take it as accurate

dsol:=dsolve({deq,ic},numeric,range=a..b);

NiM+SSVkc29sRzYiZio2I0koeF9ya2Y0NUdGJTYpSSRyZXNHRiVJJWRhdGFHRiVJJXZhcnNHRiVJKXNvbG5wcm9jR0YlSSlvdXRwb2ludEdGJUkmbmRzb2xHRiVJImlHRiU2I0lpbkNvcHlyaWdodH4oYyl+MjAwMH5ieX5XYXRlcmxvb35NYXBsZX5JbmMufkFsbH5yaWdodHN+cmVzZXJ2ZWQuR0YlRiVDKj5JOF9FbnZEU051bWVyaWNTYXZlRGlnaXRzR0YlSSdEaWdpdHNHRiU+RjYiIzpAJS9JLV9FbnZJbkZzb2x2ZUdGJUkldHJ1ZUdJKnByb3RlY3RlZEdGPT44KC0mSSZldmFsZkdGPTYjRjU2IzkkPkY/LUZCRkQ+OCVgNiRJJEdldEdGJUkkU2V0R0YlYjYjSSt0aGlzbW9kdWxlR0YlNiNJJURhdGFHRiVGJUZLRiVGJUYlRiVJK19tNTAxNjMzMjRHRiVGUT44Ji1fRklGTDYjUShzeXN2YXJzRiU+OCctRlc2I1Evc29sbl9wcm9jZWR1cmVGJUAkNC1JJXR5cGVHRj02JEY/LkkobnVtZXJpY0dGPUMkQCwtSSdtZW1iZXJHRj02JEZFNytRJnN0YXJ0RiUuSSZzdGFydEdGJVEnbWV0aG9kRiUuSSdtZXRob2RHRiVRJWxlZnRGJVEmcmlnaHRGJVEpbGVmdGRhdGFGJVEqcmlnaHRkYXRhRiVRK2VuZ2luZWRhdGFGJUMkPjgkLUZlbjYjLUkoY29udmVydEdGPTYkRkUuSSdzdHJpbmdHRj1AJi1GXG82JEZjcC5JJmFycmF5R0Y9Ty1JJWV2YWxHRj02JEZjcCIiIjBGY3BRKXByb2NuYW1lRiVPRmNwLUZjbzYkRkU3JlElbGFzdEYlLkklbGFzdEdGJVEoaW5pdGlhbEYlLkkoaW5pdGlhbEdGJUMkRmJwQCQtRlxvNiRGY3AuSSVsaXN0R0Y9TzcjLUkkc2VxR0Y9NiQvJkZVNiMsJjgqRmRxRmRxRmRxJkZjcEZecy9GYHM7IiIhRmRxMy1GXG82JEY/SSI9R0Y9LUZjbzYkLUkkbGhzR0Y9NiNGPzckRl5yRl9yQyRAJS1GXG82JC1JJHJoc0dGPUZddEZlcj5GY3AtRmVuNiMvRl5yNyRGZHMtSSNvcEdGPTYjRmN0WVFLaW5pdGlhbH52YWx1ZXN+bXVzdH5iZX5zcGVjaWZpZWR+aW5+YX5saXN0RiVGYnIvRj9RLnNvbG5wcm9jZWR1cmVGJU8tRmJxNiNGZW4vRj9GWU9GVUAlMDkhSSh1bmtub3duR0Y9Ty0uRmh1RkRDJD44KS1JKHBvaW50dG9HRj02IyYtRlc2I1Ewc29sbl9wcm9jZWR1cmVzRiU2I0Zkc08tLkZfdkZEWiVDJD5GY3AtRmVuRl10RmhyRiVZRiVGJUYlRiU=

dsol(0);

NiM3JC9JInRHNiIkIiIhRigvLUkid0dGJjYjRiUkISIiRig=

for j from 0 by 1 to steps do

print("step=",j,"TSapprox=",x[0,j],"Maple=", dsol(a+j*h));

od;

NihRJnN0ZXA9NiIiIiFRKlRTYXBwcm94PUYkISIiUSdNYXBsZT1GJDckL0kidEdGJCRGJUYlLy1JIndHRiQ2I0YrJEYnRiU=

NihRJnN0ZXA9NiIiIiJRKlRTYXBwcm94PUYkJCEwKysrKyFbdyIpISM6USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKysiRikvLUkid0dGJDYjRi0kITAiR3goMypcdyIpRik=

NihRJnN0ZXA9NiIiIiNRKlRTYXBwcm94PUYkJCEwOUtIJylIR2onISM6USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKysjRikvLUkid0dGJDYjRi0kITBSdi4qZSFIaidGKQ==

NihRJnN0ZXA9NiIiIiRRKlRTYXBwcm94PUYkJCEwPUJnb1BQSCYhIzpRJ01hcGxlPUYkNyQvSSJ0R0YkJCIwKysrKysrKyRGKS8tSSJ3R0YkNiNGLSQhMGc1JEgxI1FIJkYp

NihRJnN0ZXA9NiIiIiVRKlRTYXBwcm94PUYkJCEwRVopKUc9aDYlISM6USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKyslRikvLUkid0dGJDYjRi0kITAsKGUlZSk+O1RGKQ==

NihRJnN0ZXA9NiIiIiZRKlRTYXBwcm94PUYkJCEwT1wpbyZmKXlJISM6USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKysmRikvLUkid0dGJDYjRi0kITA8T1VjUCp5SUYp

NihRJnN0ZXA9NiIiIidRKlRTYXBwcm94PUYkJCEwOD0zemh1PCMhIzpRJ01hcGxlPUYkNyQvSSJ0R0YkJCIwKysrKysrKydGKS8tSSJ3R0YkNiNGLSQhMFZZQThSdjwjRik=

NihRJnN0ZXA9NiIiIihRKlRTYXBwcm94PUYkJCEwMCRmRWdKQDkhIzpRJ01hcGxlPUYkNyQvSSJ0R0YkJCIwKysrKysrKyhGKS8tSSJ3R0YkNiNGLSQhMGBlbHgiUkA5Rik=

NihRJnN0ZXA9NiIiIilRKlRTYXBwcm94PUYkJCEwRDwiKSlwPEMkKSEjO1EnTWFwbGU9RiQ3JC9JInRHRiQkIjArKysrKysrKSEjOi8tSSJ3R0YkNiNGLSQhMG9Ka3k0XEspRik=

NihRJnN0ZXA9NiIiIipRKlRTYXBwcm94PUYkJCEwKyRwNyV6cFclISM7USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKysqISM6Ly1JIndHRiQ2I0YtJCEwI1IqUio0b1pXRik=

NihRJnN0ZXA9NiIiIzVRKlRTYXBwcm94PUYkJCEwdXBKJEcqPi4kISM7USdNYXBsZT1GJDckL0kidEdGJCQiMCsrKysrKysiISM5Ly1JIndHRiQ2I0YtJCEwbUh4K3JFLiRGKQ==

for j from 0 by 1 to steps do

print("step=",j,"ABSerror=",abs(x[0,j]-op([2,2], dsol(a+j*h))));

od;

NiZRJnN0ZXA9NiIiIiFRKkFCU2Vycm9yPUYkJEYlRiU=

NiZRJnN0ZXA9NiIiIiJRKkFCU2Vycm9yPUYkJCIrIkd4KDM+ISM6

NiZRJnN0ZXA9NiIiIiNRKkFCU2Vycm9yPUYkJCIrRFZ1LXchIzo=

NiZRJnN0ZXA9NiIiIiRRKkFCU2Vycm9yPUYkJCIrVShHVkgpISM6

NiZRJnN0ZXA9NiIiIiVRKkFCU2Vycm9yPUYkJCIrdlJkSCEpISM6

NiZRJnN0ZXA9NiIiIiZRKkFCU2Vycm9yPUYkJCIrIm9RJip6KCEjOg==

NiZRJnN0ZXA9NiIiIidRKkFCU2Vycm9yPUYkJCIrSUc5TXghIzo=

NiZRJnN0ZXA9NiIiIihRKkFCU2Vycm9yPUYkJCIrW2wqXGQoISM6

NiZRJnN0ZXA9NiIiIilRKkFCU2Vycm9yPUYkJCIsVjkkKXpLKCEjOw==

NiZRJnN0ZXA9NiIiIipRKkFCU2Vycm9yPUYkJCIsIzRJImUsKCEjOw==

NiZRJnN0ZXA9NiIiIzVRKkFCU2Vycm9yPUYkJCIsIypmWDx5JyEjOw==

h^(methorder);