Motivating Example 

volume(r):=(4/3)*%pi*r^3;

volume(6.0);

tax(price):=price*(10/100);

tax(1000);

Motivating Example

sphere(r):=block(
area:4*%pi*r^2,
volume:(4/3)*%pi*r^3,
print("Area"=area),
print("Volume"=volume)
)$

sphere(6)$

area;

volume;

sphere(r):=block([area],
area:4*%pi*r^2,
volume:(4/3)*%pi*r^3,
print("Area"=area),
print("Volume"=volume)
)$

sphere(6)$

area;

Example 1

sol1: ode2('diff(y,x)+y=x, y, x);

solinitial:ic1(sol1, x=0, y=1);

y(x):=%e^(-x)*((x-1)*%e^x+2);

y(1/5),numer;

y(2/5),numer;

y(3/5),numer;

y(4/5),numer;

y(5/5),numer;

y

f(x,y):=x-y$
x0:0$
y0:1$
xend:1$
n:5$
h:(xend-x0)/n$
x:x0$
y:y0$
for i: 1 step 1 thru n do
block(
y:f(x,y)*h+y,
x:x+h,
print(x,y)
);

eular(f,x0,y0,xend,n):=block([h,x,y,vx,vy],
h:(xend-x0)/(n),
x:x0,
y:y0,
vx: makelist(0, n + 1),
vy: makelist(0, n + 1),
vx[1]: x0,
vy[1]: y0,
for i: 1 step 1 thru n do
block(
vy[i+1]:y:f(x,y)*h+y,
vx[i+1]:x:x+h
),
[vx, vy]
)$
[lx, ly]:eular(lambda([x, y],x-y),0,1,1,5)$

plot2d([[discrete, lx, ly],%e^(-x)*((x-1)*%e^x+2)],[x,0,1],[style, points, lines], [color, red, blue],


        [point_type, asterisk],


        [legend, "Numerical", "Exact"],[xlabel, "x"],[ylabel, "y"])$

With Step size 0.01

eular(f,x0,y0,xend,n):=block([h,x,y,vx,vy],
h:(xend-x0)/(n),
x:x0,
y:y0,
vx: makelist(0, n + 1),
vy: makelist(0, n + 1),
vx[1]: x0,
vy[1]: y0,
for i: 1 step 1 thru n do
block(
vy[i+1]:y:f(x,y)*h+y,
vx[i+1]:x:x+h
),
[vx, vy]
)$
[lx, ly]:eular(lambda([x, y],x-y),0,1,1,100)$

plot2d([[discrete, lx, ly],%e^(-x)*((x-1)*%e^x+2)],[x,0,1],[style, points, lines], [color, red, blue],


        [point_type, asterisk],


        [legend, "Numerical", "Exact"],[xlabel, "x"],[ylabel, "y"])$

Example 2

runge(f, x0, y0, x1, n) := block([h, x, y, vx, vy, k1, k2, k3, k4],
   h: bfloat((x1 - x0)/(n)),
   x: x0,
   y: y0,
   vx: makelist(0, n + 1),
   vy: makelist(0, n + 1),
   vx[1]: x0,
   vy[1]: y0,
   for i from 1 thru n do (
      k1: bfloat(h*f(x, y)),
      k2: bfloat(h*f(x + h/2, y + k1/2)),
      k3: bfloat(h*f(x + h/2, y + k2/2)),
      k4: bfloat(h*f(x + h, y + k3)),
      vy[i + 1]: y: y + (k1 + 2*k2 + 2*k3 + k4)/6,
      vx[i + 1]: x: x + h
   ),
   [vx, vy]
)$
 
[x, y]: runge(lambda([x, y], (1+y)/x), 1, 1, 10, 100)$

plot2d([discrete, x, y])$