text:
# Lab 1

## Euler's Method
### Problem 1

Solve the following IVP using Euler's Method:

$ x'(t) = -\sin(t+x)$

$ x(0) = -\pi/2$,

where the exact solution is $ x = \sin^{-1} \bigg(\frac{1-t^2}{1+t^2}\bigg) - t$

```python
'''
We are trying to solve an Ordinary Differential Equation
using the Euler's method. For this example we are taking
the example of:
dy/dt = -ky : IVP => y = y0 for t = 0
We will solve this without using the scipy
odeint package
'''

# import the basic libraries
import math
import numpy as np                # used for numerical analysis
import matplotlib.pyplot as plt   # used for plotting Graphs

# initialize the problem
x0 = -np.pi/2
t =0
x = x0
#k ,t , y =  1, 0.0, 5.0
#y0 = y
tt, yy = [], []
tf = 2                         # final value for the time

dt = 0.05                         # step size
#exact_sol = []
# describing the equation in function

def f(t,x):
    return -math.sin(t+x)

# now thw Euler method
exact_sol = []
while t<= tf:
    tt.append(t)
    yy.append(x)
    x = x+dt*f(t,x)
    exact = -math.asin((1-t**2)/(1+t**2)) - t
    exact_sol.append(exact)
    #exact = -math.asin((1-t**2)/(1+t**2))
    #exact_sol.append(exact)

    t += dt

# finding the Exact Result

#time
exact = -math.asin((1-t**2)/(1+t**2)) - t

#print(exact)
# now the Plotting part
print(len(tt))

plt.plot(tt,yy,'o', label="dt=%.4f"%(dt))
plt.plot(tt, exact_sol,'k', label="Exact Solution")
plt.xlabel("time")
plt.ylabel('y')
plt.legend(loc='best')

# plotting the absolute error

plt.legend(loc='best')

plt.title("ODE by Euler Method")
plt.show()
plt.savefig('Euler_Method.png')
```

```python
# Euler's Method

f = lambda t,x : -np.sin(t + x)
n = 100
k = 2/n
t = np.linspace(0,2,n)
x = np.zeros(n)
x[0] = -np.pi/2

for i in range(n-1):
  x[i+1] = x[i] + k*f(t[i],x[i])

plt.plot(t,x,'b')
plt.plot(t,-np.arcsin((1-t*t)/(1+t*t))-t,'g')
```

```python
# Python Code to find approximation
# of a ordinary differential equation
# using euler method.

# Consider a differential equation
xx = []
yy = []
def func( t, y ):
    return (-np.sin(x+y))

# Function for euler formula
def euler( x0, y, h, x ):
    temp = -0

    # Iterating till the point at which we
    # need approximation
    while x0 < x:
        xx.append(x)
        yy.append(y)
        temp = y
        y = y + h * func(x0, y)
        x0 = x0 + h

    # Printing approximation
    print("Approximate solution at time = ", x, " is ", "%.6f"% y)

# Driver Code
# Initial Values
x0 = 0
y0 = -np.pi/2
h = 0.025

# Value of x at which we need approximation
x = 2

euler(x0, y0, h, x)
```

    Approximate solution at time =  2  is  -1.943906

## Taylor-Series Method of order 4

### Problem 2

Solve the following IVP in the interval $[-1,1]$ using Taylor-Series method of order $4$:

$ x'  = \cos(t) - \sin(x) + t^2$

$ x(-1) = 3$

```python
# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them

# we will begin by defining the function
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return math.cos(t) - math.sin(x) + t**2

# now we differentiate the same function and find the
# next orders of derivatives
def f2(t,x):
    return -math.sin(t) - math.cos(x)*f1(t,x) + 2*t

def f3(t,x):
    return -math.cos(t) + 2 - f2(t,x) * math.cos(x) + math.sin(x)*f1(t,x)**2

def f4(t,x):
    return math.sin(t) - (f2(t,x)*f1(t,x)*(-math.sin(x)) + math.cos(x)* f3(t,x)) + (math.cos(x)* f1(t,x)**3 + math.sin(x) * 2* f2(t,x) )

# also the function that finds the Error

def error(dt,m):
    return ((dt**m)/factorial(m+1)) * (xx[len(xx)-1] - xx[len(xx)-2])

# and the exact solution to the problem(for reference)
#def exact(t):
#    return math.sqrt(np.pi)*math.erf(t)*0.5

# Now we write down the values known to us
m = 4                       # order parameter
dt = 0.09                    # step Size
t,x = -1 ,3                 # initial values of both
tf = 1                      # final point to reach
#xe = 0                      # for the error
# initialised arrays for graph plotting
tt = []
xx = []
#xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    #xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x) + np.power(dt,4)/24*f4(t,x)       # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1

# Plotting and Visualization

#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'o',linewidth=2,color='red', label ='Taylor Method Order = 4')
#plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method.png')

```

## Taylor Series method of order 3

### Problem 3

Solve the following delay-differential equation using the Taylor-Series method of order 3:
$ x'(t) = x(t-1)$; $(t\ge0)$

$ x(t) = t^2 $ ; $(-1 \le t \le 0)$

```python
# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them

# we will begin by defining the function
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return (t-1)**2

# now we differentiate the same function and find the
# next orders of derivatives
def f2(t,x):
    return 2*(t-1)

def f3(t,x):
    return 2

# and the exact solution to the problem(for reference)
#def exact(t):
#    return math.sqrt(np.pi)*math.erf(t)*0.5

# Now we write down the values known to us
m = 3                       # order parameter
dt = 0.05                    # step Size
t,x = 0 ,0                 # initial values of both
tf = 2                      # final point to reach
#xe = 0                      # for the error
# initialised arrays for graph plotting
tt = []
xx = []
#xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    #xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x)        # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1

# Plotting and Visualization

#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'--',linewidth=2,color='red', label ='Taylor Method Order = 3')
#plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method_order_3.png')

```

## Taylor-Series Method of order 4
### Problem 4

Use the Taylor-Series method of order $4$ to compute $\int e^{-t^2} ds$ by solving the IVP on the interval $ t \in [0,2]$

$ x' = e^{-t^2}$

$x(0) = 0 $

```python
'''
Use the Taylor Series method of order 4 to compute $\int_0^2 e^-s^2 ds
by solving the Initial Value problem on the interval t \in [0,2]
Given dx/dt = e^-t^2
and   x(0) = 0
Check from the table of error function x(2) = (pi^(1/2)/2). erf(2) = 0.8820813907
We need to find the derivative till the fourth order of the Given function, and add them to a Taylor
Series to form the solution
'''

# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them

# we will begin by defining the function
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return exp((-t**2))

# now we differentiate the same function and find the
# next orders of derivatives
def f2(t,x):
    return -2*t*f1(t,x)

def f3(t,x):
    return 4*np.power(t,2)*f1(t,x) - 2*f1(t,x)

def f4(t,x):
    return 12*t*f1(t,x) - 8*np.power(t,3)*f1(t,x)

# also the function that finds the Error

def error(dt,m):
    return ((dt**m)/factorial(m+1)) * (xx[len(xx)-1] - xx[len(xx)-2])

# and the exact solution to the problem(for reference)
def exact(t):
    return math.sqrt(np.pi)*math.erf(t)*0.5

# Now we write down the values known to us
m = 4                       # order parameter
dt = 0.04                    # step Size
t,x = 0 ,0                    # initial values of both
tf = 2                      # final point to reach
xe = 0
# initialised arrays for graph plotting
tt = []
xx = []
xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x) + np.power(dt,4)/24*f4(t,x)       # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1

# Plotting and Visualization

import scipy

from scipy.special import erf
t = np.linspace(0,2,500)
#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'o',linewidth=2,color='red', label ='Taylor Method Order = 4')
plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method3.png')

# checking for the value at t = 2
# writing all the values

print("The Value that x(t) takes for different t:\n", xx)

# printing the error

print("The Error Value is:\n ",error(dt,m))

#  difference between exact value and Taylor approximated

print("The Difference between the Exact solution at t = 2 and Taylor approximated is:\n", math.erf(2)-xx[len(xx)-1])

```

## Runge-Kutta method of Order 4 && Rungw-Kutta Fehlberg Method

### Problem 5
 Solve the following IVP on the interval $[1,3]$

 $ x' = t^{-2} (tx - x^2)$

 $x(1) = 2$

 using the Runge Kutta Method and the adaptive runge kutta method

Solve the IVP on the interval $[1,3]$

$ x' = t^{-2} (tx - x^2)$

$ x(1) = 2$

```python
import math
import matplotlib.pyplot as plt

# function
dy = lambda t,x: (t**(-2))*(t*x - x**2)
f = lambda t: ((0.5 + ln(t))**(-1))*t

#initialization

ti = 1
tf = 3
h = 0.1
n = int((tf-ti)/h)
t = 1
x = 2

# visualization
print('t  \t\t x \t\t f(t)')
print('%f \t %f \t %f'% (t,x,f(t)))

t_plot = []
x_RK4 = []
x_analytical = []

# RK4 method
for i in range(1, n+1):
  t_plot.append(t)
  x_RK4.append(x)
  x_analytical.append(f(t))

  # Euler method
  k1 = dy(t,x)
  k2 = dy(t+h/2, x + k1*h/2)
  k3 = dy(t+h/2, x + k2*h/2)
  k4 = dy(t+h, x+ k3*h)

  # new estimation
  x = x + 1/6 *(k1 + 2*k2 + 2*k3 +k4) *h
  # increase by the step size
  t = t+h
  print('%f \t %f \t %f'% (t,x,f(t)))

# visualization
t_plot.append(t)
x_RK4.append(x)
x_analytical.append(f(t))

plt.plot(t_plot, x_RK4, 'ro', t_plot, x_analytical)
plt.xlabel('x')
plt.ylabel('y')
plt.legend(["RK4","Analytical"])

```

```python
# RK2/4
import numpy as np
import matplotlib.pyplot as plt
n = 256
t = np.linspace(1,3,256)
h = 2/n
x = np.zeros(n)
x[0] = 2
f = lambda t,x: (t*x - x*x)/(t*t)

for i in range(n-1):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h,x[i] + h*k1)
  x[i+1] = x[i] + h/2*(k1+k2)#RK2
#for i in range(n-1):
#  k1 = f(t[i],x[i])
 # k2 = f(t[i]+h/2,x[i] + h*k1/2)
 # k3 = f(t[i]+h/2,x[i] + h*k2/2)
 # k4 = f(t[i]+h,x[i] + h*k3)
#  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)#RK4

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t)))
```

## Runge-Kutta-Fehlberg (RKF45) method

```python
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define min(a, b) ((a) < (b) ? (a) : (b))

double f (double t, double y);

int main ()
{
	double x0,y0,k1,k2,k3,k4,k5,k6,w,w1,w2,y,x,xn,y1,y2,h,R,epsilon,delta, g;
	int i=0;

	//boundary conditions
	epsilon = 0.00001;
	x0=1;
	y0=2;
	xn=3;
	h=0.2;
	x=x0;
	y=y0;

	printf("Step %d: x=%f\t y=%f\n", i, x, y);
	while (x<=xn)
	{
		h = min(h, 3-x);

		k1 = h*f(x,y);
		k2 = h*f(x+h/4, y+k1/4);
		k3 = h*f(x+3*h/8,  y+3*k1/32+9*k2/32);
		k4 = h*f(x+12*h/13, y+1932*k1/2197-7200*k2/2197+7296*k3/2197);
		k5 = h*f(x+h, y+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
		k6 = h*f(x+h/2, y-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);

		y1 = y + (25*k1/216+1408*k3/2565+2197*k4/4104-k5/5);
		y2 = y + (16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55);

		R = ((fabs(y2-y1)))/h;
		delta = 0.84*(pow((epsilon/R), 0.25));

		if (R<=epsilon)
		{
			x = x+h;
			y = y1;
			i = i+1;
			printf("Step %d: x=%f\t y=%f\n", i, x, y);
			h = delta*h;
		}
		else
		{
			h = delta*h;
            x = x+h;
		}
	}
}

double f(double x, double y)
{
	double v;
	v = pow(x,-2)*(x*y - pow(y,2));
	return v;
}
```

```python
% Runge-Kutta-Fehlberg method
 % Approximate the solution to the initial-value problem
 %  ; 1<=t<=3 ; y(0)=0.5; tolerance=10^(-5); hmin=0.01; hmax=0.25

 f = @(t,y) (t^(-2)*(t*y - y^2));
 a = input('Enter left end ponit, a:  ');
 b = input('Enter right end point, b:  ');
 n = input('Enter no. of subintervals, n: ');
 alpha = input('Enter the initial condition, alpha:  ');
 tol = input('Enter tolerance, tol:  ');
 hmin = input('Enter minimum mesh, hmin:  ');
 hmax = input('Enter maximum mesh, hmax:  ');

 h = hmax;
 t = a;
 w = alpha;
 fprintf('  t         w\n');
 fprintf('%5.4f  %11.8f  \n', t, w);

 while t < b
    k1 = h*f(t,w);
    k2 = h*f(t+h/4,w+k1/4);k3 = h*f(t+3*h/8,w+(3*k1+9*k2)/32);
    k4 = h*f(t+12*h/13,w+(1932*k1-7200*k2+7296*k3)/2197);
    k5 = h*f(t+h,w+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
    k6 = h*f(t+h/2,w-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
    r = abs(k1/360-128*k3/4275-2197*k4/75240.0+k5/50+2*k6/55)/h;
    if r <= tol
       t = t+h;
       w = w+25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
       fprintf('%5.4f  %11.8f  \n', t, w);
       plot(t,w,'b*'); grid on;
       xlabel('t values'); ylabel('w values');
       hold on;
    end

    if r > 1.0e-20
       delta = 0.84 * exp(0.25 * log(tol / r));
    else
       delta = 10.0;
    end

    if delta <= 0.1
       h = 0.1*h;
    else
       if delta >= 4
          h = 4.0 * h;
       else
          h = delta * h;
       end
    end

   if h > hmax
      h = hmax;
   end
   if h < hmin
     fprintf('Minimal h exceeded\n');
   else
      if t+h > b
         if abs(b-t) < tol
            t = b;
         else
            h = b-t;
         end
      end
   end
 end
```

# Lab 2

Solving via all the methods

## ABM4 Method

Solve the following IVP using ABM4 method in the interval [1,3]

 Solve the following IVP on the interval $[1,3]$

 $ x' = t^{-2} (tx - x^2)$

```python
from numpy import *
import numpy as np
from math import *
from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

'''
Write down the function that we need to solve for
x' = t^{-2}(tx-x^2)
'''

def f(t,x):
    return (t**(-2))*(t*x - (x**2))

'''
another derivative of the function function
'''

def dfx(t,x):
    return (1/np.power(t,2))*(t-2*x)

# analytical solution

# backward Euler
def beuler(t0, tn, n, x0):
    #h = abs(tn - t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0

    for k in range(0,n):
        err = 1
        zold = x[k] + h*f(t[k], x[k])
        I = 0

        while err > 10**(-10) and I < 5:
            F = x[k] + h*f(t[k+1], zold) - zold
            dF = h*dfx(t[k+1], zold)-1
            znew = zold - F/dF
            err = abs(znew- zold)
            I += 1
        x[k+1] = znew
    return x

def RK4(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0
    for i in range(0,n):
        K1 = f(t[i], x[i])
        K2 = f(t[i]+h/2, x[i]+K1*h/2)
        K3 = f(t[i] + h/2, x[i]+K2*h/2)
        K4 = f(t[i] + h, x[i]+K3*h)

        x[i+1] = x[i] + h*(K1+2*K2+ 2*K3+K4)/6
    return x

def AdBash3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    return x

## ABM4 pred corr

def PreCorr3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i] + h*(9*K0 + 19*K1 - 5*K2 + K3)/24

    return x

fg = 1
n= 200
t0 = 1
h = 0.015
tn = 3
x0 = 2
t = linspace(t0, tn , n+1)
xe = AdBash3(t0,tn, n, x0)
xb = beuler(t0,tn,n,x0)
xpc = PreCorr3(t0, tn, n,x0)

# print
#exact = sol(t0, tn ,n , x0)

plot(t,xe,'o',color='yellow',label ='Adam Bashforth 3')
plot(t, xb, '--', color='red', label ='Backward Euler')
plot(t,xpc,'cyan',label='Predictor Corrector')
#plot(t,(t* (1/2 + log(t))**(-1)),'blue',label ='Exact Solution')

#t = linspace(t0, tn, 401)
#xsol = sol(t,t0,x0)

#plot(t,xsol, color='green', label='Exact')

##title('n = %d' %n)
t = t0
exact = []
time = []
while(t<=tn):

    m = (t* (1/2 + log(t))**(-1))

    t = t+0.006666666666666821
    exact.append(m)
    time.append(t)

plt.plot(time,exact,'blue',label='Exact')
#axis([0,tn, -60,40])
legend(loc='upper right')
plt.savefig('ADM_h=0.015.png')
```

```python
#ABM4
f = lambda t,x: (t*x - x*x)/(t*t)
n = 100
h = 2/n
t = np.linspace(1,3,n)
x = np.zeros(n)
x[0] = 2

for i in range(3):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h/2,x[i] + h*k1/2)
  k3 = f(t[i]+h/2,x[i] + h*k2/2)
  k4 = f(t[i]+h,x[i] + h*k3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)

for i in range(3,n-1):
  fn = f(t[i],x[i])
  fn1 = f(t[i-1],x[i-1])
  fn2 = f(t[i-2],x[i-2])
  fn3 = f(t[i-3],x[i-3])
  p = x[i] + h/24*(55*fn -59*fn1 + 37*fn2 -9*fn3)
  fnp = f(t[i],p)
  x[i+1] = x[i] + h/24*(9*fnp + 19*fn -5*fn1 + fn2)

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t)))
```

```python
from numpy import *
from matplotlib.pyplot import *
from matplotlib.patches import Patch
#f in the IVP y’ = f(t,y), y(t0)=y0
def f(t,y):
  return (t-3.2)*y + 8*t*exp((t-3.2)**2/2)*cos(4*t**2)

def dfy(t,y):
  return t-3.2
#analytic solution to the IVP y’ = f(t,y), y(t0)=y0

def sol(t,t0,y0):
  C = y0*exp(-(t0-3.2)**2/2)-sin(4*t0**2)
  return exp((t-3.2)**2/2)*(sin(4*t**2) + C)

#Backwards Euler (Implicit) Method
def beuler(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(0,tn,n+1)
  y = zeros(n+1)
  y[0] = y0
  for k in range(0,n):
    err = 1
    zold = y[k] + h*f(t[k],y[k]) #Use Forward Euler for initial guess
    I = 0
    #Use Newton’s Method to solve implicit equation for y[k+1]
    while err > 10**(-10) and I < 5: #NM is limited to 5 iterations
      F = y[k] + h*f(t[k+1],zold)-zold
      dF = h*dfy(t[k+1],zold)-1
      znew = zold - F/dF
      err = abs(znew-zold)
      zold = znew
      I +=1
    y[k+1] = znew
  return y

#Runge-Kutta "Classic" Order 4 method
def RK4(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  y[0] = y0
  for i in range(0,n):
    K1 = f(t[i],y[i])
    K2 = f(t[i]+h/2,y[i]+K1*h/2)
    K3 = f(t[i]+h/2,y[i]+K2*h/2)
    K4 = f(t[i]+h,y[i]+K3*h)
    y[i+1] = y[i] + h*(K1+2*K2+2*K3+K4)/6
  return y

#Adams-Bashforth 3 Step Method
def AdBash3(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  y[0:3] = RK4(t0,t0+2*h,2,y0)
  K1 = f(t[1],y[1])
  K2 = f(t[0],y[0])
  for i in range(2,n):
    K3 = K2
    K2 = K1
    K1 = f(t[i],y[i])
    y[i+1] = y[i] + h*(23*K1-16*K2+5*K3)/12
  return y

#Adams-Bashforth 3/Moulton 4 Step Predictor/Corrector
def PreCorr3(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  #Calculate initial steps with Runge-Kutta 4
  y[0:3] = RK4(t0,t0+2*h,2,y0)
  K1 = f(t[1],y[1])
  K2 = f(t[0],y[0])
  for i in range(2,n):
    K3 = K2
    K2 = K1
    K1 = f(t[i],y[i])
    #Adams-Bashforth Predictor
    y[i+1] = y[i] + h*(23*K1-16*K2+5*K3)/12
    K0 = f(t[i+1],y[i+1])
    #Adams-Moulton Corrector
    y[i+1] = y[i] + h*(9*K0+19*K1-5*K2+K3)/24
  return y

#Script to produce graphs
fg =1
n = 300
t0 = 0
tn = 6
y0 = .75
t = linspace(t0,tn,n+1)
ye = AdBash3(t0,tn,n,y0)
yb = beuler(t0,tn,n,y0)
ypc = PreCorr3(t0,tn,n,y0)
figure(fg)
plot(t,ye,'red',label='Adams-Bashforth 3')
plot(t,yb,'black',label='Backward Euler')
plot(t,yb,'cyan',label='Predictor/Corrector 3/4')
t = linspace(t0,tn,401)
ysol = sol(t,t0,y0)
plot(t,ysol,color ='green',label='Exact')
title('n = %d' % n)
axis([0,tn,-60,40])
legend(loc='lower left')

```

## Milne Simpson Method

```python
from numpy import *
import numpy as np
from math import *
from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

'''
Write down the function that we need to solve for
x' = t^{-2}(tx-x^2)
'''

def f(t,x):
    return (t**(-2))*(t*x - (x**2))

'''
another derivative of the function function
'''

def dfx(t,x):
    return (1/np.power(t,2))*(t-2*x)

# analytical solution

# backward Euler
def beuler(t0, tn, n, x0):
    #h = abs(tn - t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0

    for k in range(0,n):
        err = 1
        zold = x[k] + h*f(t[k], x[k])
        I = 0

        while err > 10**(-10) and I < 5:
            F = x[k] + h*f(t[k+1], zold) - zold
            dF = h*dfx(t[k+1], zold)-1
            znew = zold - F/dF
            err = abs(znew- zold)
            I += 1
        x[k+1] = znew
    return x

def RK4(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0
    for i in range(0,n):
        K1 = f(t[i], x[i])
        K2 = f(t[i]+h/2, x[i]+K1*h/2)
        K3 = f(t[i] + h/2, x[i]+K2*h/2)
        K4 = f(t[i] + h, x[i]+K3*h)

        x[i+1] = x[i] + h*(K1+2*K2+ 2*K3+K4)/6
    return x

def AdBash3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    return x

## ABM4 pred corr

def PreCorr3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i] + h*(9*K0 + 19*K1 - 5*K2 + K3)/24

    return x

def Milne_Simspon(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i-3] + 4*h*(2*K3 - K2 + 2*K1)/3
        #x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i-1] + h*(K0 + 4*K3 + K2)/3

    return x

fg = 1
n= 200
t0 = 1
h = 0.015
tn = 3
x0 = 2
t = linspace(t0, tn , n+1)
xe = AdBash3(t0,tn, n, x0)
xb = beuler(t0,tn,n,x0)
xpc = PreCorr3(t0, tn, n,x0)
xxpc = Milne_Simspon(t0, tn, n, x0)

# print
#exact = sol(t0, tn ,n , x0)

plot(t,xe,'o',color='yellow',label ='Adam Bashforth 3')
plot(t, xb, '--', color='red', label ='Backward Euler')
plot(t,xpc,'cyan',label='Predictor Corrector')
plot(t,xxpc,'red',label='Milne Simspon')
#plot(t,(t* (1/2 + log(t))**(-1)),'blue',label ='Exact Solution')

#t = linspace(t0, tn, 401)
#xsol = sol(t,t0,x0)

#plot(t,xsol, color='green', label='Exact')

##title('n = %d' %n)
t = t0
exact = []
time = []
while(t<=tn):

    m = (t* (1/2 + log(t))**(-1))

    t = t+0.006666666666666821
    exact.append(m)
    time.append(t)

plt.plot(time,exact,'blue',label='Exact')
#axis([0,tn, -60,40])
legend(loc='upper right')
#plt.savefig('ADM_h=0.015.png')
```

```python
#Milne-Simpson

f = lambda t,x: (t*x - x*x)/(t*t)
n = 100
h = 2/n
t = np.linspace(1,3,n)
x = np.zeros(n)
x[0] = 2

for i in range(3):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h/2,x[i] + h*k1/2)
  k3 = f(t[i]+h/2,x[i] + h*k2/2)
  k4 = f(t[i]+h,x[i] + h*k3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)

for i in range(3,n-1):
  fn = f(t[i],x[i])
  fn1 = f(t[i-1],x[i-1])
  fn2 = f(t[i-2],x[i-2])
  p = x[i-3] + 4*h/3*(2*fn -fn1 + 2*fn2)
  fnp = f(t[i],p)
  x[i+1] = x[i-1] + h/3*(fnp + 4*fn + fn1)

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t)))
```

## System of ODEs using Runge-Kutta 4

```python
import math
import time                         # import time to use for performance analysis
import numpy as np                  # import numpy for array space
import matplotlib.pyplot as plt     # import matplotlib for graphing functions
from scipy.integrate import odeint  # import scipy to use the ordinary differencial equation integral function

# Ordinary Differential Equation
def dy_dx(y, x):                    # takes inputs y and x
    return x + 2*y      # returns the function y/(e^x - 1)

# Runge-Kutta Forumla based off the formula used in class
def rungeKutta(y, x, h):                                # takes y, x, and h inputs
    k1 = dy_dx(y, x)                                    # solves k1 using the differential equation function
    k2 = dy_dx(y + ((0.5 * h) * k1), x + 0.5 * h)       # solves k2 based on the answer from k1 using the differential equation function
    k3 = dy_dx(y + ((0.5 * h) * k2), x + 0.5 * h)       # solves k3 based on the answer from k2 using the differential equation function
    k4 = dy_dx(y + (h * k3), x + h)                     # solves k4 based on the answer from k3 using the differential equation function

    t4 = (1.0 / 6.0) * (k1 + (2 * k2) + (2 * k3) + k4)  # solves for t4 by taking a 6th of k1 + 2*k2 + 2*k3 + k4

    y = y + (t4 * h)                                    # solves for y by taking the initial y value and adding it to t4*h
    return y                                            # returns y

# initial variable values
y = 6
           # initial y value of 5 as outlined in the assignment
x = 0           # initial x value of 1 as outlined in the assignment
h = 0.02        # initial h value of 0.02 as outlined in the assignment
n = 2000        # initial n value of 2000 chosen between 1000 or 2000 as outlined in the assignment

# Runge-Kutta x and y arrays
xxsr = []        # x-space runge-kutta array to store the values of x for the runge-kutta function
yysr = []        # y-space runge-kutta array to store the values of y for the runge-kutta function

# ODEint x and y arrays, solution, and time analysis
tso = time.time()                               # time start for ODEint function solution

xso = np.linspace(0.0, 0.2 , 30)       # x-space ODEint from 1 to n*h (40) plus 1 with a step size of n
yyso = odeint(dy_dx, y, xso)                     # y-space ODEint useing the odeint function from scicpy to find the y-space

teo = time.time()                               # time end fore ODEint function solution
tto = teo - tso                                 # total time the ODEint function to solve

# graphing ODEint
plt.title("ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'g-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Exact Solution.png')                                                  # displays the graph

# Runge-Kutta solution and time analysis
tsr = time.time()                               # time start for runge-kutta function solution

while (x <= 0.2):                           # for loop to run the runge-kutta function n number of times
    xxsr.append(x)                               # append the x value to the x-space runge-kutta array
    yysr.append(y)                               # append the y value to the y-space runge-kutta array

    y = rungeKutta(y, x, h)                     # update the y value using the rungeKutta function

    x += h                                      # update the x value by moving one step forward (0.02)

ter = time.time()                               # time end for runge-kutta function solution
ttr = ter - tsr                                 # total time the runge-kutta function to solve

td = ttr - tto                                  # time difference between ODEint and runge-kutta function

# graphing runge-kutta
plt.title("Runge-Kutta Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xxsr, yysr, 'b--', label = "Runge Kutta")             # set the runge-kutta to be blue and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Runge_Kutta.png')                                                  # displays the graph

# solutions
print("\nODEint Solution:            ", yyso[-1])    # ODEint function solution
print("Runge-Kutta Solution:        ", yysr[-1])     # Runge-Kutta function solution

# Print statement for time difference
print("\nODEint Time:         ", tto)               # print the ODEint time
print("Runge Kutta Time:    ", ttr)                 # print the runge-kutta time
print("ODEint is ", td, " seconds faster\n\n")      # print the difference between ODEint and runge-kutta

# error calculation
error = 0                                                           # initial error value of 0
errorRange = []                                                     # array to store error over xn
errorSpace = np.linspace(0.0,0.2, 30)                   # error space for error analysis
for i in range(len(yysr)):                                   # for loop to run through every x values
    error += (np.abs(yysr[i] - yyso[i])/yyso[i]) * 100                 # sum all the error values using the percentage error formula
    errorRange.append((np.abs(yysr[i] - yyso[i])/yyso[i]) * 100)
    print("Percent Error at x =", i, ":", (np.abs(yysr[i] - yyso[i])/yyso[i]) * 100)   # print error at each x value

print("\nAverage Error Percent:", error/((int)(n * h) + 1), "\n")     # print the total error divided by the total number of x values

# graphing error
#plt.title("Error Analysis")                                     # set the title of the graph
#plt.xlabel("xn")                                                # set the x label on the graph
#plt.ylabel("error")                                             # set the y label on the graph
#plt.plot(errorSpace, errorRange, label = "Error over Xn")       # create the line and label it
#plt.legend()                                                    # shows the legend on the graph
#plt.show()                                                      # displays the graph
#
# graphing both functions
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xxsr, yysr, 'bo', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison.png')
```

```python
import math
import time                         # import time to use for performance analysis
import numpy as np                  # import numpy for array space
import matplotlib.pyplot as plt     # import matplotlib for graphing functions
from scipy.integrate import odeint  # import scipy to use the ordinary differencial equation integral function

# Ordinary Differential Equation
def dy_dx(y, x):                    # takes inputs y and x
    return 3*x + 2*y      # returns the function y/(e^x - 1)

# Runge-Kutta Forumla based off the formula used in class
def rungeKutta(y, x, h):                                # takes y, x, and h inputs
    k1 = dy_dx(y, x)                                    # solves k1 using the differential equation function
    k2 = dy_dx(y + ((0.5 * h) * k1), x + 0.5 * h)       # solves k2 based on the answer from k1 using the differential equation function
    k3 = dy_dx(y + ((0.5 * h) * k2), x + 0.5 * h)       # solves k3 based on the answer from k2 using the differential equation function
    k4 = dy_dx(y + (h * k3), x + h)                     # solves k4 based on the answer from k3 using the differential equation function

    t4 = (1.0 / 6.0) * (k1 + (2 * k2) + (2 * k3) + k4)  # solves for t4 by taking a 6th of k1 + 2*k2 + 2*k3 + k4

    y = y + (t4 * h)                                    # solves for y by taking the initial y value and adding it to t4*h
    return y                                            # returns y

# initial variable values
y = 4
           # initial y value of 5 as outlined in the assignment
x = 0           # initial x value of 1 as outlined in the assignment
h = 0.02        # initial h value of 0.02 as outlined in the assignment
n = 2000        # initial n value of 2000 chosen between 1000 or 2000 as outlined in the assignment

# Runge-Kutta x and y arrays
xsr = []        # x-space runge-kutta array to store the values of x for the runge-kutta function
ysr = []        # y-space runge-kutta array to store the values of y for the runge-kutta function

# ODEint x and y arrays, solution, and time analysis
tso = time.time()                               # time start for ODEint function solution

xso = np.linspace(0.0, 0.2 , 30)       # x-space ODEint from 1 to n*h (40) plus 1 with a step size of n
yso = odeint(dy_dx, y, xso)                     # y-space ODEint useing the odeint function from scicpy to find the y-space

teo = time.time()                               # time end fore ODEint function solution
tto = teo - tso                                 # total time the ODEint function to solve

# graphing ODEint
plt.title("ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yso, 'g-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Exact Solution1.png')                                                  # displays the graph

# Runge-Kutta solution and time analysis
tsr = time.time()                               # time start for runge-kutta function solution

while (x <= 0.2):                           # for loop to run the runge-kutta function n number of times
    xsr.append(x)                               # append the x value to the x-space runge-kutta array
    ysr.append(y)                               # append the y value to the y-space runge-kutta array

    y = rungeKutta(y, x, h)                     # update the y value using the rungeKutta function

    x += h                                      # update the x value by moving one step forward (0.02)

ter = time.time()                               # time end for runge-kutta function solution
ttr = ter - tsr                                 # total time the runge-kutta function to solve

td = ttr - tto                                  # time difference between ODEint and runge-kutta function

# graphing runge-kutta
plt.title("Runge-Kutta Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xsr, ysr, 'b--', label = "Runge Kutta")             # set the runge-kutta to be blue and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Runge_Kutta1.png')                                                  # displays the graph

# solutions
print("\nODEint Solution:            ", yso[-1])    # ODEint function solution
print("Runge-Kutta Solution:        ", ysr[-1])     # Runge-Kutta function solution

# Print statement for time difference
print("\nODEint Time:         ", tto)               # print the ODEint time
print("Runge Kutta Time:    ", ttr)                 # print the runge-kutta time
print("ODEint is ", td, " seconds faster\n\n")      # print the difference between ODEint and runge-kutta

# error calculation
error = 0                                                           # initial error value of 0
errorRange = []                                                     # array to store error over xn
errorSpace = np.linspace(0.0,0.2, 30)                   # error space for error analysis
for i in range(len(ysr)):                                   # for loop to run through every x values
    error += (np.abs(ysr[i] - yso[i])/yso[i]) * 100                 # sum all the error values using the percentage error formula
    errorRange.append((np.abs(ysr[i] - yso[i])/yso[i]) * 100)
    print("Percent Error at x =", i, ":", (np.abs(ysr[i] - yso[i])/yso[i]) * 100)   # print error at each x value

print("\nAverage Error Percent:", error/((int)(n * h) + 1), "\n")     # print the total error divided by the total number of x values

# graphing error
#plt.title("Error Analysis")                                     # set the title of the graph
#plt.xlabel("xn")                                                # set the x label on the graph
#plt.ylabel("error")                                             # set the y label on the graph
#plt.plot(errorSpace, errorRange, label = "Error over Xn")       # create the line and label it
#plt.legend()                                                    # shows the legend on the graph
#plt.show()                                                      # displays the graph
#
# graphing both functions
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xsr, ysr, 'bo', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison1.png')
```

```python
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xso, yso, 'b-', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison1.png')
```

```python
#RK4 for a system of equations

fx = lambda t,x,y: x + 2*y
fy = lambda t,x,y: 3*x + 2*y
n = 100
h = 0.2/n
t = np.linspace(0,0.2,n)
x = np.zeros(n)
y = np.zeros(n)
x[0] = 6
y[0] = 4

for i in range(n-1):
  k1 = fx(t[i],x[i],y[i])
  m1 = fy(t[i],x[i],y[i])
  k2 = fx(t[i]+h/2,x[i] + h*k1/2,y[i] + h*m1/2)
  m2 = fy(t[i]+h/2,x[i] + h*k1/2,y[i] + h*m1/2)
  k3 = fx(t[i]+h/2,x[i] + h*k2/2,y[i] + h*m2/2)
  m3 = fy(t[i]+h/2,x[i] + h*k2/2,y[i] + h*m2/2)
  k4 = fx(t[i]+h,x[i] + h*k3,y[i] + h*m3)
  m4 = fy(t[i]+h,x[i] + h*k3,y[i] + h*m3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)
  y[i+1] = y[i] + h/6*(m1+2*m2+2*m3+m4)

plt.plot(t,x)
plt.plot(t,4*np.exp(4*t)+2*np.exp(-t))
plt.plot(t,y)
plt.plot(t,6*np.exp(4*t)-2*np.exp(-t))
```

```python
import matplotlib.pyplot as plt

x = 6
y = 4
t = 0

h = 0.05

while (t < 1):
  m1 = (x + 2*y)
  k1 = 3*x + 2*y
  ft2 = t + (h/2)
  fx2 = x+ (h/2) * m1
  fy2 = y + (h/2) * k1
  m2 = (2 * fx2) + (4 * fy2)

  gt2 = t + (h/2)
  gx2 = x + (h/2) * m1
  gy2 = y + (h/2) * k1
  k2 = -(gx2)+ (6 *gy2)

  ft3 = (t + (h/2))
  fx3 = x +(h/2)* m2
  fy3 = y + (h/2)*k2
  m3 = (2 * fx3) + (4* fy3)

  gt3 = t + (h/2)
  gx3 = x + (h/2) * m2
  gy3 = y + (h/2) * k2
  k3 = -(gx3)+ (6 *gy3)

  ft4 = (t + h)
  fx4 = x +h* m3
  fy4 = y + h*k3
  m4 = (2 * fx4) + (4* fy4)

  gt4 = t + h
  gx4 = x + h * m3
  gy4 = y + h * k3
  k4 = -(gx4)+ (6 *gy4)

  t = t+h

  x = x +( h/6) * (m1 + (2 * m2) + (2 * m3) + m4)
  y = y +( h/6) * (k1 + (2 * k2) + (2 * k3) + k4)

  print(t, '\t', '\t', x, '\t', y)
  plt.plot(t,x,'rs')
  plt.plot(t,y,'bs')

plt.show()

```

```python
import numpy as np
from math import sqrt
import matplotlib.pyplot as plt
import numba
import time
start_time = time.clock()

a=1
b=1
c=1
d=1

# Equations:
@numba.jit()

#du/dt=V(u,t)
def V(u,t):
    x, y, vx, vy = u
    return np.array([vy,vx,a*x+b*y,c*x+d*y])

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = t[1]-t[0]
    for i in range(n):
        k1 = h * f(u[i], t[i])
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
    return u, t

u, t  = rk4(V, np.array([1., 0., 1. , 0.]) , 0. , 10. , 100000)
x,y, vx,vy  = u.T
# plt.plot(t, x, t,y)
plt.semilogy(t, x, t,y)
plt.grid('on')
print("Execution time:",time.clock() - start_time, "seconds")
plt.show()
```

```python

```

```python

```

## Linear Shooting Method

```python
#Linear Shooting Method

fy = lambda t,x,y: 2*t/(1+t*t)*y - 2/(1+t*t)*x + 1
fx = lambda t,x,y: y
n = 100
t = np.linspace(0,4,n)
h = (t[n-1]-t[0])/n
x = np.zeros(n)
u = np.zeros(n)
du = np.zeros(n)
v = np.zeros(n)
dv = np.zeros(n)
u[0] = 1.25
dv[0] = 1

for i in range(n-1):
  k1 = fx(t[i],u[i],du[i])
  m1 = fy(t[i],u[i],du[i])
  k2 = fx(t[i]+h/2,u[i] + h*k1/2,du[i] + h*m1/2)
  m2 = fy(t[i]+h/2,u[i] + h*k1/2,du[i] + h*m1/2)
  k3 = fx(t[i]+h/2,u[i] + h*k2/2,du[i] + h*m2/2)
  m3 = fy(t[i]+h/2,u[i] + h*k2/2,du[i] + h*m2/2)
  k4 = fx(t[i]+h,u[i] + h*k3,du[i] + h*m3)
  m4 = fy(t[i]+h,u[i] + h*k3,du[i] + h*m3)
  u[i+1] = u[i] + h/6*(k1+2*k2+2*k3+k4)
  du[i+1] = du[i] + h/6*(m1+2*m2+2*m3+m4)

for i in range(n-1):
  k1 = fx(t[i],v[i],dv[i])
  m1 = fy(t[i],v[i],dv[i])
  k2 = fx(t[i]+h/2,v[i] + h*k1/2,dv[i] + h*m1/2)
  m2 = fy(t[i]+h/2,v[i] + h*k1/2,dv[i] + h*m1/2)
  k3 = fx(t[i]+h/2,v[i] + h*k2/2,dv[i] + h*m2/2)
  m3 = fy(t[i]+h/2,v[i] + h*k2/2,dv[i] + h*m2/2)
  k4 = fx(t[i]+h,v[i] + h*k3,dv[i] + h*m3)
  m4 = fy(t[i]+h,v[i] + h*k3,dv[i] + h*m3)
  v[i+1] = v[i] + h/6*(k1+2*k2+2*k3+k4)
  dv[i+1] = dv[i] + h/6*(m1+2*m2+2*m3+m4)

w = (-0.95 - u[n-1])/(v[n-1])*v
x = u + w

plt.plot(t,x)
```

```python
# matlab

 p = @(x) ((2*x)/(1+x^2));
 q = @(x) (-2/1 + x^2);
 r = @(x) (1);

 a=0; b=4; alpha =1.25; beta = -0.95;  n = 50;

 fprintf('   x          w\n');
 h = (b-a)/n;
 u1 = alpha;
 u2 = 0;
 v1 = 0;
 v2 = 1;
 u = zeros(2,n);
 v = zeros(2,n);

 for i = 1 : n
   x = a+(i-1)*h;
   t = x+0.5*h;
   k11 = h*u2;
   k12 = h*(p(x)*u2+q(x)*u1+r(x));
   k21 = h*(u2+0.5*k12);
   k22 = h*(p(t)*(u2+0.5*k12)+q(t)*(u1+0.5*k11)+r(t));
   k31 = h*(u2+0.5*k22);
   k32 = h*(p(t)*(u2+0.5*k22)+q(t)*(u1+0.5*k21)+r(t));
   t = x+h;
   k41 = h*(u2+k32);
   k42 = h*(p(t)*(u2+k32)+q(t)*(u1+k31)+r(t));
   u1 = u1+(k11+2*(k21+k31)+k41)/6;
   u2 = u2+(k12+2*(k22+k32)+k42)/6;
   k11 = h*v2;
   k12 = h*(p(x)*v2+q(x)*v1);
   t = x+0.5*h;
   k21 = h*(v2+0.5*k12);
   k22 = h*(p(t)*(v2+0.5*k12)+q(t)*(v1+0.5*k11));
   k31 = h*(v2+0.5*k22);
   k32 = h*(p(t)*(v2+0.5*k22)+q(t)*(v1+0.5*k21));
   t = x+h;
   k41 = h*(v2+k32);
   k42 = h*(p(t)*(v2+k32)+q(t)*(v1+k31));
   v1 = v1+(k11+2*(k21+k31)+k41)/6;
   v2 = v2+(k12+2*(k22+k32)+k42)/6;
   u(1,i) = u1;
   u(2,i) = u2;
   v(1,i) = v1;
   v(2,i) = v2;
 end

 w1 = alpha;
 z = (beta-u(1,n))/v(1,n);
 x = a;
 i = 0;
 fprintf('%5.4f   %11.8f \n', x, w1);

 for i = 1 : n
   x = a+i*h;
   w1 = u(1,i)+z*v(1,i);
   w2 = u(2,i)+z*v(2,i);
   fprintf('%5.4f   %11.8f \n', x, w1);
   plot(x,w1,'b*'); grid on;
   xlabel('t values'); ylabel('w values');
   hold on;
 end
```

```python
import numpy

#-----------------------------------------------------------------------------

def shoot( f, a, b, z1, z2, t, tol ):
    """Implements the shooting method to solve second order BVPs
    USAGE:
        y = shoot(f, a, b, z1, z2, t, tol)
    INPUT:
        f     - function dy/dt = f(y,t).  Since we are solving a second-
                order boundary-value problem that has been transformed
                into a first order system, this function should return a
                1x2 array with the first entry equal to y and the second
                entry equal to y'.
        a     - solution value at the left boundary: a = y(t[0]).
        b     - solution value at the right boundary: b = y(t[n-1]).
        z1    - first initial estimate of y'(t[0]).
        z1    - second initial estimate of y'(t[0]).
        t     - array of n time values to determine y at.
        tol   - allowable tolerance on right boundary: | b - y[n-1] | < tol
    OUTPUT:
        y     - array of solution function values corresponding to the
                values in the supplied array t.
    NOTE:
        This function assumes that the second order BVP has been converted to
        a first order system of two equations.  The secant method is used to
        refine the initial values of y' used for the initial value problems.
    """

    from diffeq import rk4

    max_iter = 25   # Maximum number of shooting iterations

    n = len( t )    # Determine the size of the arrays we will generate

    # Compute solution to first initial value problem (IVP) with y'(a) = z1.
    # Because we are using the secant method to refine our estimates of z =
    # y', we don't really need all the solution of the IVP, just the last
    # point of it -- this is saved in w1.

    y = rk4( f, [a,z1], t )
    w1 = y[n-1,0]

    print("%2d: z = %10.3e, error = %10.3e" % ( 0, z1, b - w1 ))

    # Begin the main loop.  We will compute the solution of a second IVP and
    # then use the both solutions to refine our estimate of y'(a).  This
    # second solution then replaces the first and a new "second" solution is
    # generated.  This process continues until we either solve the problem to
    # within the specified tolerance or we exceed the maximum number of
    # allowable iterations.

    for i in range( max_iter ):

        # Solve second initial value problem, using y'(a) = z2.  We need to
        # retain the entire solution vector y since if y(t(n)) is close enough
        # to b for us to stop then the first column of y becomes our solution
        # vector.

        y = rk4( f, [a,z2], t )
        w2 = y[n-1,0]

        print( "%2d: z = %10.3e, error = %10.3e" % ( i+1, z2, b - w2 ))

        # Check to see if we are done...

        if abs( b - w2 ) < tol:
            break

        # Compute the new approximations to the initial value of the first
        # derivative.  We compute z2 using a linear fit through (z1,w1) and
        # (z2,w2) where w1 and w2 are the estimates at t=b of the initial
        # value problems solved above with y1'(a) = z1 and y2'(a) = z2.  The
        # new value for z1 is the old value of z2.

        #z1, z2 = ( z2, z1 + ( z2 - z1 ) / ( w2 - w1 ) * ( b - w1 ) )
        z1, z2 = ( z2, z2 + ( z2 - z1 ) / ( w2 - w1 ) * ( b - w2 ) )
        w1 = w2

    # All done.  Check to see if we really solved the problem, and then return
    # the solution.

    if abs( b - w2 ) >= tol:
        print ("\a**** ERROR ****")
        print ("Maximum number of iterations (%d) exceeded" % max_iter)
        print ("Returned values may not have desired accuracy")
        print ("Error estimate of returned solution is %e" % ( b - w2 ))

    return y[:,0]

#-----------------------------------------------------------------------------

def fd( u, v, w, t, a, b ):
    """Implements the shooting method to solve linear second order BVPs
    Compute finite difference solution to the BVP
        x'' = u(t) + v(t) x + w(t) x'
        x(t[0]) = a, x(t[n-1]) = b
    t should be passed in as an n element array.   u, v, and w should be
    either n element arrays corresponding to u(t), v(t) and w(t) or
    scalars, in which case an n element array with the given value is
    generated for each of them.
    USAGE:
        x = fd(u, v, w, t, a, b)
    INPUT:
        u,v,w - arrays containing u(t), v(t), and w(t) values.  May be
                specified as Python lists, NumPy arrays, or scalars.  In
                each case they are converted to NumPy arrays.
        t     - array of n time values to determine x at
        a     - solution value at the left boundary: a = x(t[0])
        b     - solution value at the right boundary: b = x(t[n-1])
    OUTPUT:
        x     - array of solution function values corresponding to the
                values in the supplied array t.
    """

    # Get the dimension of t and make sure that t is an n-element vector

    if type( t ) != numpy.ndarray:
        if type( t ) == list:
            t = numpy.array( t )
        else:
            t = numpy.array( [ float( t ) ] )

    n = len( t )

    # Make sure that u, v, and w are either scalars or n-element vectors.
    # If they are scalars then we create vectors with the scalar value in
    # each position.

    if type( u ) == int or type( u ) == float:
        u = numpy.array( [ float( u ) ] * n )

    if type( v ) == int or type( v ) == float:
        v = numpy.array( [ float( v ) ] * n )

    if type( w ) == int or type( w ) == float:
        w = numpy.array( [ float( w ) ] * n )

    # Compute the stepsize.  It is assumed that all elements in t are
    # equally spaced.

    h = t[1] - t[0];

    # Construct tridiagonal system; boundary conditions appear as first and
    # last equations in system.

    A = -( 1.0 + w[1:n] * h / 2.0 )
    A[-1] = 0.0

    C = -( 1.0 - w[0:n-1] * h / 2.0 )
    C[0] = 0.0

    D = 2.0 + h * h * v
    D[0] = D[n-1] = 1.0

    B = - h * h * u
    B[0] = a
    B[n-1] = b

    # Solve tridiagonal system

    for i in range( 1, n ):
        xmult = A[i-1] / D[i-1]
        D[i] = D[i] - xmult * C[i-1]
        B[i] = B[i] - xmult * B[i-1]

    x = numpy.zeros( n )
    x[n-1] = B[n-1] / D[n-1]

    for i in range( n - 2, -1, -1 ):
        x[i] = ( B[i] - C[i] * x[i+1] ) / D[i]

    return x

#-----------------------------------------------------------------------------

if __name__ == "__main__":

    import math
    from pylab import *

    # Solves x'' = x + 4exp(t), x(0)=1, x(1/2) = 2exp(1/2) using both the
    # finite difference method and the shooting method.

    # Set up interval.  We will solve the problem for both n=64 and n=128.

    a = 1.25
    b = -0.95
    n1 = 64
    n2 = 128
    t1 = linspace( a, b, n1 )
    t2 = linspace( a, b, n2 )

    # Compute exact solutions.  The transpose of the solution is taken because
    # both the finite difference function fd() and the shooting method function
    # shoot() return nx1 vectors rather than 1xn vectors.

    def exact(t):
        return exp(t) * ( 1 + 2 * t )

    x1 = exact( t1 )
    x2 = exact( t2 )

    # Compute finite difference solutions

    xfd1 = fd( 4 * exp( t1 ), 1, 0, t1, 1, 2 * exp( 0.5 ) )
    xfd2 = fd( 4 * exp( t2 ), 1, 0, t2, 1, 2 * exp( 0.5 ) )

    # Compute shooting method solutions

    def f(x,t):
        return array( [x[1], x[0]+4*exp(t)] )

    xs1 = shoot( f, exp(a), 2*exp(b), 3.0, 4.0, t1, 1e-5 )
    xs2 = shoot( f, exp(a), 2*exp(b), 3.0, 4.0, t2, 1e-5 )

    # Prepare for display; set interactive mode to true so each plot
    # is shown as it is generated

    interactive( True )

    # Plot solutions

    plot( t1, xfd1, 'ro', t2, xfd2, 'b-' )
    title( 'Finite Difference Method' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='lower right' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla()
    plot( t1, xs1, 'ro', t2, xs2, 'b-' )
    title( 'Shooting Method' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='lower right' )
    draw()
    z = input( "Press ENTER to continue..." )

    # Plot errors

    cla()
    plot( t1, xfd1 - x1, 'ro', t2, xfd2 - x2, 'b-' )
    title( 'Finite Difference Errors' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='center' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla()
    plot( t1, xs1 - x1, 'ro', t2, xs2 - x2, 'b-' )
    title( 'Shooting Method Errors' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='center' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla()
    plot( t1, xfd1 - x1, 'ro-', t1, xs1 - x1, 'b-' )
    # fiddle with ymin value so that shooting method error is not at
    # bottom of window
    win = list( axis() )
    win[2] = win[2] - 0.1 * ( win[3] - win[2] );
    axis( tuple( win ) )
    title( 'Finite Difference and Shooting Method Errors: %d points' % n1 )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( 'Finite Differences', 'Shooting' ), loc='center' )
    draw()
    z = input( "Press ENTER to quit..." );
```

```python
# diffeq
from __future__ import print_function
import numpy as np

#-----------------------------------------------------------------------------

def euler( f, x0, t ):
    """Euler's method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = euler(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """

    n = len( t )
    x = np.array( [x0] * n )
    for i in range( n - 1 ):
        x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )

    return x

#-----------------------------------------------------------------------------

def heun( f, x0, t ):
    """Heun's method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = heun(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """
    n = len( t )
    x = np.array( [x0] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + k1, t[i+1] )
        x[i+1] = x[i] + ( k1 + k2 ) / 2.0

    return x

#-----------------------------------------------------------------------------

def rk2a( f, x0, t ):
    """Second-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk2a(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Analysis", 6th Edition, by Burden and Faires, Brooks-Cole, 1997.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] ) / 2.0
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )

    return x

#-----------------------------------------------------------------------------

def rk2b( f, x0, t ):
    """Second-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk2b(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Mathematics and Computing" 4th Edition, by Cheney and Kincaid,
        Brooks-Cole, 1999.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + k1, t[i+1] )
        x[i+1] = x[i] + ( k1 + k2 ) / 2.0

    return x

#-----------------------------------------------------------------------------

def rk4( f, x0, t ):
    """Fourth-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk4(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + 0.5 * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + k3, t[i+1] )
        x[i+1] = x[i] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0

    return x

#-----------------------------------------------------------------------------

def rk45( f, x0, t ):
    """Fourth-order Runge-Kutta method with error estimate.
    USAGE:
        x, err = rk45(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
        err   - NumPy array containing estimate of errors at each step.  If
                a system is being solved, err will be an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Mathematics and Computing" 6th Edition, by Cheney and Kincaid,
        Brooks-Cole, 2008.
    """

    # Coefficients used to compute the independent variable argument of f

    c20  =   2.500000000000000e-01  #  1/4
    c30  =   3.750000000000000e-01  #  3/8
    c40  =   9.230769230769231e-01  #  12/13
    c50  =   1.000000000000000e+00  #  1
    c60  =   5.000000000000000e-01  #  1/2

    # Coefficients used to compute the dependent variable argument of f

    c21 =   2.500000000000000e-01  #  1/4
    c31 =   9.375000000000000e-02  #  3/32
    c32 =   2.812500000000000e-01  #  9/32
    c41 =   8.793809740555303e-01  #  1932/2197
    c42 =  -3.277196176604461e+00  # -7200/2197
    c43 =   3.320892125625853e+00  #  7296/2197
    c51 =   2.032407407407407e+00  #  439/216
    c52 =  -8.000000000000000e+00  # -8
    c53 =   7.173489278752436e+00  #  3680/513
    c54 =  -2.058966861598441e-01  # -845/4104
    c61 =  -2.962962962962963e-01  # -8/27
    c62 =   2.000000000000000e+00  #  2
    c63 =  -1.381676413255361e+00  # -3544/2565
    c64 =   4.529727095516569e-01  #  1859/4104
    c65 =  -2.750000000000000e-01  # -11/40

    # Coefficients used to compute 4th order RK estimate

    a1  =   1.157407407407407e-01  #  25/216
    a2  =   0.000000000000000e-00  #  0
    a3  =   5.489278752436647e-01  #  1408/2565
    a4  =   5.353313840155945e-01  #  2197/4104
    a5  =  -2.000000000000000e-01  # -1/5

    b1  =   1.185185185185185e-01  #  16.0/135.0
    b2  =   0.000000000000000e-00  #  0
    b3  =   5.189863547758284e-01  #  6656.0/12825.0
    b4  =   5.061314903420167e-01  #  28561.0/56430.0
    b5  =  -1.800000000000000e-01  # -9.0/50.0
    b6  =   3.636363636363636e-02  #  2.0/55.0

    n = len( t )
    x = np.array( [ x0 ] * n )
    e = np.array( [ 0 * x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + c21 * k1, t[i] + c20 * h )
        k3 = h * f( x[i] + c31 * k1 + c32 * k2, t[i] + c30 * h )
        k4 = h * f( x[i] + c41 * k1 + c42 * k2 + c43 * k3, t[i] + c40 * h )
        k5 = h * f( x[i] + c51 * k1 + c52 * k2 + c53 * k3 + c54 * k4, \
                        t[i] + h )
        k6 = h * f( \
            x[i] + c61 * k1 + c62 * k2 + c63 * k3 + c64 * k4 + c65 * k5, \
            t[i] + c60 * h )

        x[i+1] = x[i] + a1 * k1 + a3 * k3 + a4 * k4 + a5 * k5
        x5 = x[i] + b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5 + b6 * k6

        e[i+1] = abs( x5 - x[i+1] )

    return ( x, e )

#-----------------------------------------------------------------------------

def rkf( f, a, b, x0, tol, hmax, hmin ):
    """Runge-Kutta-Fehlberg method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        t, x = rkf(f, a, b, x0, tol, hmax, hmin)
    INPUT:
        f     - function equal to dx/dt = f(x,t)
        a     - left-hand endpoint of interval (initial condition is here)
        b     - right-hand endpoint of interval
        x0    - initial x value: x0 = x(a)
        tol   - maximum value of local truncation error estimate
        hmax  - maximum step size
        hmin  - minimum step size
    OUTPUT:
        t     - NumPy array of independent variable values
        x     - NumPy array of corresponding solution function values
    NOTES:
        This function implements 4th-5th order Runge-Kutta-Fehlberg Method
        to solve the initial value problem
           dx
           -- = f(x,t),     x(a) = x0
           dt
        on the interval [a,b].
        Based on pseudocode presented in "Numerical Analysis", 6th Edition,
        by Burden and Faires, Brooks-Cole, 1997.
    """

    # Coefficients used to compute the independent variable argument of f

    a2  =   2.500000000000000e-01  #  1/4
    a3  =   3.750000000000000e-01  #  3/8
    a4  =   9.230769230769231e-01  #  12/13
    a5  =   1.000000000000000e+00  #  1
    a6  =   5.000000000000000e-01  #  1/2

    # Coefficients used to compute the dependent variable argument of f

    b21 =   2.500000000000000e-01  #  1/4
    b31 =   9.375000000000000e-02  #  3/32
    b32 =   2.812500000000000e-01  #  9/32
    b41 =   8.793809740555303e-01  #  1932/2197
    b42 =  -3.277196176604461e+00  # -7200/2197
    b43 =   3.320892125625853e+00  #  7296/2197
    b51 =   2.032407407407407e+00  #  439/216
    b52 =  -8.000000000000000e+00  # -8
    b53 =   7.173489278752436e+00  #  3680/513
    b54 =  -2.058966861598441e-01  # -845/4104
    b61 =  -2.962962962962963e-01  # -8/27
    b62 =   2.000000000000000e+00  #  2
    b63 =  -1.381676413255361e+00  # -3544/2565
    b64 =   4.529727095516569e-01  #  1859/4104
    b65 =  -2.750000000000000e-01  # -11/40

    # Coefficients used to compute local truncation error estimate.  These
    # come from subtracting a 4th order RK estimate from a 5th order RK
    # estimate.

    r1  =   2.777777777777778e-03  #  1/360
    r3  =  -2.994152046783626e-02  # -128/4275
    r4  =  -2.919989367357789e-02  # -2197/75240
    r5  =   2.000000000000000e-02  #  1/50
    r6  =   3.636363636363636e-02  #  2/55

    # Coefficients used to compute 4th order RK estimate

    c1  =   1.157407407407407e-01  #  25/216
    c3  =   5.489278752436647e-01  #  1408/2565
    c4  =   5.353313840155945e-01  #  2197/4104
    c5  =  -2.000000000000000e-01  # -1/5

    # Set t and x according to initial condition and assume that h starts
    # with a value that is as large as possible.

    t = a
    x = x0
    h = hmax

    # Initialize arrays that will be returned

    T = np.array( [t] )
    X = np.array( [x] )

    while t < b:

        # Adjust step size when we get to last interval

        if t + h > b:
            h = b - t;

        # Compute values needed to compute truncation error estimate and
        # the 4th order RK estimate.

        k1 = h * f( x, t )
        k2 = h * f( x + b21 * k1, t + a2 * h )
        k3 = h * f( x + b31 * k1 + b32 * k2, t + a3 * h )
        k4 = h * f( x + b41 * k1 + b42 * k2 + b43 * k3, t + a4 * h )
        k5 = h * f( x + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4, t + a5 * h )
        k6 = h * f( x + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5, \
                    t + a6 * h )

        # Compute the estimate of the local truncation error.  If it's small
        # enough then we accept this step and save the 4th order estimate.

        r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
        if len( np.shape( r ) ) > 0:
            r = max( r )
        if r <= tol:
            t = t + h
            x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
            T = np.append( T, t )
            X = np.append( X, [x], 0 )

        # Now compute next step size, and make sure that it is not too big or
        # too small.

        h = h * min( max( 0.84 * ( tol / r )**0.25, 0.1 ), 4.0 )

        if h > hmax:
            h = hmax
        elif h < hmin:
            print( "Error: stepsize should be smaller than %e." % hmin )
            break

    # endwhile

    return ( T, X )

#-----------------------------------------------------------------------------

def pc4( f, x0, t ):
    """Adams-Bashforth-Moulton 4th order predictor-corrector method
    USAGE:
        x = pc4(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This function used the Adams-Bashforth-Moulton predictor-corrector
        method to solve the initial value problem
            dx
            -- = f(x,t),     x(t(1)) = x0
            dt
        at the t values stored in the t array (so the interval of solution is
        [t[0], t[N-1]].  The 4th-order Runge-Kutta method is used to generate
        the first three values of the solution.  Notice that it works equally
        well for scalar functions f(x,t) (in the case of a single 1st order
        ODE) or for vector functions f(x,t) (in the case of multiple 1st order
        ODEs).
    """

    n = len( t )
    x = np.array( [ x0 ] * n )

    # Start up with 4th order Runge-Kutta (single-step method).  The extra
    # code involving f0, f1, f2, and f3 helps us get ready for the multi-step
    # method to follow in order to minimize the number of function evaluations
    # needed.

    f1 = f2 = f3 = 0
    for i in range( min( 3, n - 1 ) ):
        h = t[i+1] - t[i]
        f0 = f( x[i], t[i] )
        k1 = h * f0
        k2 = h * f( x[i] + 0.5 * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + k3, t[i+1] )
        x[i+1] = x[i] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0
        f1, f2, f3 = ( f0, f1, f2 )

    # Begin Adams-Bashforth-Moulton steps

    for i in range( 3, n - 1 ):
        h = t[i+1] - t[i]
        f0 = f( x[i], t[i] )
        w = x[i] + h * ( 55.0 * f0 - 59.0 * f1 + 37.0 * f2 - 9.0 * f3 ) / 24.0
        fw = f( w, t[i+1] )
        x[i+1] = x[i] + h * ( 9.0 * fw + 19.0 * f0 - 5.0 * f1 + f2 ) / 24.0
        f1, f2, f3 = ( f0, f1, f2 )

    return x

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

if __name__ == "__main__":
    from pylab import *

    def f( x, t ):
        return x * np.sin( t )

    a, b = ( 0.0, 10.0 )
    x0 = -1.0

    n = 51
    t = np.linspace( a, b, n )
    h = t[1] - t[0];
    tol = 1e-6

    # compute various numerical solutions
    x_euler = euler( f, x0, t )
    x_heun = heun( f, x0, t )
    x_rk2 = rk2a( f, x0, t )
    x_rk4 = rk4( f, x0, t )
    x_pc4 = pc4( f, x0, t )
    t_rkf, x_rkf = rkf( f, a, b, x0, tol, 1.0, 0.01 ) # unequally spaced t

    # compute true solution values in equal spaced and unequally spaced cases
    x = -np.exp( 1.0 - np.cos( t ) )
    xrkf = -np.exp( 1.0 - np.cos( t_rkf ) )

#   figure( 1 )
    subplot( 2, 2, 1 )
    plot( t, x_euler, 'b-o', t, x_heun, 'g-o', t, x_rk2, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    title( 'Solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)' % h )
    legend( ( 'Euler $O(h)$', 'Heun $O(h^2)$', 'Runge-Kutta $O(h^2)$' ),
            loc='lower left' )

#   figure( 2 )
    subplot( 2, 2, 2 )
    plot( t, x_euler - x, 'b-o', t, x_heun - x, 'g-o', t, x_rk2 - x, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x - x^*$' )
    title( 'Errors in solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)'
           % h )
    legend( ( 'Euler $O(h)$', 'Heun $O(h^2)$', 'Runge-Kutta $O(h^2)$' ),
            loc='upper left' )

#   figure( 3 )
    subplot( 2, 2, 3 )
    plot( t, x_rk4, 'b-o', t, x_pc4, 'g-o', t_rkf, x_rkf, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    title( 'Solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)' % h )
    legend( ( 'Runge-Kutta $O(h^4)$', 'Predictor-Corrector $O(h^4)$', \
              'Runge-Kutta-Fehlberg' ), loc='lower left' )

#   figure( 4 )
    subplot( 2, 2, 4 )
    plot( t, x_rk4 - x, 'b-o', t, x_pc4 - x, 'g-o', t_rkf, x_rkf - xrkf, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x - x^*$' )
    title( 'Errors in solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)'
           % h )
    legend( ( 'Runge-Kutta $O(h^4)$', 'Predictor-Corrector $O(h^4)$', \
              'Runge-Kutta-Fehlberg' ), loc='lower left' )

    show()
```

#### Finite Differences

```python
#Finite Differences Method

n = 50
t = np.linspace(0,4,n)
h = 4/n
p = 2*t/(1+t*t)
q = -2/(1+t*t)
A = np.zeros((n-1,n-1))
B = np.zeros(n-1)
x = np.zeros(n)

for i in range(n-2):
  A[i,i] = 2 + h*h*q[i]
  A[i,i+1] = h/2 *p[i] - 1
  A[i+1,i] = -h/2*p[i+1] -1
  B[i] = -h*h
A[n-2,n-2] = 2 + h*h*q[n-2]
B[0] = -h*h + (h/2*p[0] +1)*1.25
B[n-2] = -h*h + (-h/2*p[n-2]+1)*(-0.95)
x = np.linalg.solve(A,B)

plt.plot(t[:-1],x)
```

# Lab 3

### Problem 1

Use Finite Difference method to solve the wave equation
$ u_{tt}(x,t) = 4 u_{xx}(x,t)$
for

$0 < x <1 $

```python
# Finite Difference Wave Eqaution
from mpl_toolkits.mplot3d import Axes3D
h = 0.1
k = 0.05
c = 2
r = c*k/h
n = 10
u = np.zeros((n,n))
t = np.arange(0,0.5,k)
x = np.arange(0,1,h)
A = np.zeros((n-2,n-2))
B = np.zeros(n-2)
u[:,0] = np.sin(np.pi*x) + np.sin(2*np.pi*x)
u[:,1] = u[:,0]
for i in range(n-3):
  A[i,i] = 0
  A[i,i+1] = 1
  A[i+1,i] = 1
#for i in range(1,n-1):
 # u[i,1] = (u[i+1,0] + u[i-1,0])/2
#for i in range(1,n-1):
#  u[1:n-1,i+1] = np.dot(A,u[1:n-1,i]) - u[1:n-1,i-1]
for j in range(1,n-1):
  for i in range(1,n-1):
    u[i,j+1] = u[i-1,j] + u[i+1,j] -u[i,j-1]
X,Y = np.meshgrid(x,t)
Z = np.sin(np.pi*X)*np.cos(2*np.pi*Y) + np.sin(2*np.pi*X)*np.cos(4*np.pi*Y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z,'-o',linewidth=0.75)
#ax.view_init(45, 120)

```

```python
clc;
clear;
close all;

nx = 10;
nt = 30;
a = 0;
b = 1;
t0 = 0;
tf = 0.5;

dx = (b-a)/(nx-1);
dt = (tf- t0)/(nt-1);

x = a:dx:b;
t = t0:dt:tf;

s = dt^2/dx^2;
c = 2;

%% AnalyticalSolution

UA = zeros(nx,nt);

for i = 1:nx
    for j = 1:nt
        UA(i,j) = sin(pi*x(i)) * cos(2*pi*t(j)) + sin(2*pi*x(i)) * cos(4*pi*t(j));
    end
end

figure()
contourf(UA,200,'linecolor', 'non')
title('Analytical Solution')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar

%% Numerical Solution

UN = zeros(nx, nt);

UN(:,1) = sin(pi*x) + sin(2*pi*x);
UN(:,2) = (sin(pi*x) + sin(2*pi*x)) * (1+dt);

for j = 2:nt-1
    for i = 2:nx-1
        UN (i, j+1) = (4)*s*(UN(i-1,j) - 2*UN(i,j) + UN(i+1,j)) + 2*UN(i,j) - UN(i, j-1);
    end
end

figure()
contourf(UN,200,'linecolor', 'non')
title(' Numerical Solution')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar

%% Error
E = abs(UA - UN);

figure()
contourf(E,200,'linecolor', 'non')
title(' Error')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar

```

### Problem 2

```python
from numpy import *
import numpy as np, matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display

alpha = 0.5
dx = 0.0050251
dt = 0.0050251
rho = 0.01
ten = 40.
c = sqrt(ten/rho)     # Density, tension
c1 = c
ratio =  0.5           # CFL criterium = 1
xi = np.zeros((101,3), float)                  # Declaration
k = range(0,101)

def Initialize():                            # Initial conditions
   for i in range(0, 61):
       xi[i, 0] = 0.00125*i
   for i in range (61, 101):
       xi[i, 0] = 0.1 - 0.005*(i - 80)

def animate(num):
    for i in range(1, 100):
       xi[i,2] = ratio*(xi[i+1,1]+ xi[i-1,1] - 2*xi[i,1]) + alpha*dt**2/dx*(xi[i+1,1] - xi[i,1]) - xi[i,0] + 2*xi[i,1]
    line.set_data(k,xi[k,2])          # Data to plot ,x,y
    for m in range (0,101):
       xi[m, 0] = xi[m, 1]            # Recycle array
       xi[m, 1] = xi[m, 2]
    return line

Initialize()                         # Plot initial string
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 101),
	ylim=(-15, 15))
ax.grid()                                     # Plot  grid
plt.title("Vibrating String")
line, = ax.plot(k, xi[k,0], lw=2)
for i in range(1,100):
    xi[i,1] = ratio*(xi[i+1,0]+ xi[i-1,0] - 2*xi[i,0]) + alpha*dt**2/dx*(xi[i+1,0] - xi[i,0]) - xi[i,0]
ani = animation.FuncAnimation(fig, animate,1)    # Dummy 1

writervideo = animation.FFMpegWriter(fps=60)
ani.save('wave.mp4', writer=writervideo)
plt.show()
print("finished")
```

```python
#FTCS

D = 1
h = 0.2
k = 0.02
r = 2*D*k/(h**2)
t = np.linspace(0,0.2,10)
x = np.linspace(0,1,5)
A = np.zeros((4,4))

```

```python
#Laplace

n = 5
x = np.linspace(0,4,5)
y = np.linspace(0,4,5)
a = np.diagflat(-4*np.ones(3)) + np.diagflat(np.ones(2),1) + np.diagflat(np.ones(2),-1)
A = np.kron(np.eye(3),a) + np.diagflat(np.ones(6),3) + np.diagflat(np.ones(6),-3)
u = [-100,-20,-20,-80,0,0,-260,-180,-180]
p = np.linalg.solve(A,u)
p
```

    array([ 55.71428571,  43.21428571,  27.14285714,  79.64285714,
            70.        ,  45.35714286, 112.85714286, 111.78571429,
            84.28571429])

```python
# Laplace Neumann

n = 5
x = np.linspace(0,4,5)
y = np.linspace(0,4,5)
a = np.diagflat(-4*np.ones(3)) + np.diagflat(np.ones(2),1) + np.diagflat(np.ones(2),-1)
A = np.kron(np.eye(4),a) + np.diagflat([2,2,2,1,1,1,1,1,1],3) + np.diagflat(np.ones(9),-3)
u = [-80,0,0,-80,0,0,-80,0,0,-260,-180,-180]
p = np.linalg.solve(A,u)
p
```

    array([ 71.82182413,  56.85429572,  32.23419526,  75.2165004 ,
            61.68058174,  36.04124267,  87.36359573,  78.61028818,
            50.25019366, 115.62759433, 115.14678158,  86.34924381])

```python
# SOR
from mpl_toolkits.mplot3d import Axes3D
#import plotly.express as px
n = 10
x = np.linspace(0,4,n)
y = np.linspace(0,4,n)
u = 70*np.ones((n,n))
w = 1.44646
r = 1

for i in range(1,n-2):
  u[i,n-1] = 180
  u[i,0] = 20
  u[n-1,i]=0
  u[0,i] = 80
u[1,1] = 50
u[9,1] = 10
u[1,9] = 130
u[9,9] = 90
count = 0
while np.absolute(r)>0.001:
  for i in range(1,n-2):
    for j in range(1,n-2):
      r = (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j])/4
      u[i,j] = u[i,j] + w*r
      count = count + 1

X,Y = np.meshgrid(x,y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,u,'-o',linewidth=0.75)
#ax.view_init(-140,35)
```

```python
# SOR Neumann
from mpl_toolkits.mplot3d import Axes3D
#import plotly.express as px
n = 10
x = np.linspace(0,4,n)
y = np.linspace(0,4,n)
u = 70*np.ones((n,n))
w = 1.44646
r = 1

for i in range(1,n-2):
  u[i,n-1] = 180
  u[n-1,i]=0
  u[0,i] = 80
u[1,1] = 50
u[9,1] = 10
u[1,9] = 130
u[9,9] = 90
count = 0
while np.absolute(r)>0.001:
  for i in range(1,n-2):
    r = (u[i+1,0] + u[i-1,0] + 2*u[i,1] - 4*u[i,0])/4
    for j in range(1,n-2):
      r = (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j])/4
      u[i,j] = u[i,j] + w*r
      count = count + 1

X,Y = np.meshgrid(x,y)
#fig = plt.figure()
#ax = Axes3D(fig)
#ax.plot_wireframe(X,Y,u,'-o',linewidth=0.75)
plt.pcolormesh(x,y,u)
#ax.view_init(-140,35)
```

```python

```

# Lab 4

#### Heat Equation solve

using MATLAB

```python
clc
clear
close all

nx = 10;
nt = 10;
a = 0;
b = 1;

t0 = 0;
tf = 0.2;

dx = (b-a)/(nx-1);
dt = (tf- t0)/(nt-1);

x = a:dx:b;
t = t0:dt:tf;

s = dt/dx^2;

%% Numerical
UN = zeros (nx, nt);

% Initial

UN(:,1) = 4*x - 4*(x.*x);

for j = 1:nt-1
    for i = 2:nx-1
        UN(i,j+1) = s*UN (i-1,j) + (1-2*s)*UN(i,j) + s*UN(i+1,j);
    end
end
figure()
contourf(UN,200,'linecolor', 'non')
title('Numerical Solution')
xlabel('x')
ylabel('t')
colormap(jet(256))
colorbar
caxis([0,1])

```

## DFT and FFT

```python
import numpy as np

def dft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

"""
We can ensure our implementation is correct by comparing
the results with those obtained from numpy’s fft function.
"""
x = np.random.random(1024)
np.allclose(dft(x), np.fft.fft(x))

# check time
%timeit dft(x)
%timeit np.fft.fft(x)

```

    10 loops, best of 5: 125 ms per loop
    The slowest run took 4.98 times longer than the fastest. This could mean that an intermediate result is being cached.
    10000 loops, best of 5: 16.1 µs per loop

```python
def fft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if N % 2 > 0:
        raise ValueError("must be a power of 2")
    elif N <= 2:
        return dft(x)
    else:
        X_even = fft(x[::2])
        X_odd = fft(x[1::2])
        terms = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + terms[:int(N/2)] * X_odd,
                               X_even + terms[int(N/2):] * X_odd])

x = np.random.random(1024)
np.allclose(fft(x), np.fft.fft(x))

%timeit dft(x)
%timeit fft(x)
%timeit np.fft.fft(x)
```

    10 loops, best of 5: 125 ms per loop
    10 loops, best of 5: 14.2 ms per loop
    The slowest run took 6.52 times longer than the fastest. This could mean that an intermediate result is being cached.
    100000 loops, best of 5: 16.1 µs per loop

We define another function to compute the Fourier Transform. Only this time around, we make use of vector operations instead of recursion

```python
def fft_v(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if np.log2(N) % 1 > 0:
        raise ValueError("must be a power of 2")

    N_min = min(N, 2)

    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))
    while X.shape[0] < N:
        X_even = X[:, :int(X.shape[1] / 2)]
        X_odd = X[:, int(X.shape[1] / 2):]
        terms = np.exp(-1j * np.pi * np.arange(X.shape[0])
                        / X.shape[0])[:, None]
        X = np.vstack([X_even + terms * X_odd,
                       X_even - terms * X_odd])
    return X.ravel()

```

we can ensure we obtained the correct results by comparing them with those from the numpy library.

```python
x = np.random.random(1024)
np.allclose(fft_v(x), np.fft.fft(x))
```

    True

```python
%timeit fft(x)
%timeit fft_v(x)
%timeit np.fft.fft(x)
```

    100 loops, best of 5: 14.5 ms per loop
    The slowest run took 4.33 times longer than the fastest. This could mean that an intermediate result is being cached.
    1000 loops, best of 5: 371 µs per loop
    The slowest run took 9.66 times longer than the fastest. This could mean that an intermediate result is being cached.
    100000 loops, best of 5: 16.2 µs per loop

## Forward Time Difference Heat Equation

```python
## matlab

L = 0.1;  % wall thickness
n = 10;   % number of simulation nodes
T0 = 0;  % initial temperature of the wall
T1s = 40;  % surface 1 temp
T2s = 20;  % surface 2 temperature

dx = L/n;   % node thickness

alpha = 0.0001;   % thermal diffusivity

t_final = 60 ;   % simulation time

dt = 0.1; % fixed time step

x = dx/2:dx:L-dx/2;

T = ones(n,1)*T0;

dTdt = zeros(n,1);

t = 0:dt:t_final;

for j = 1:length(t)
    for i = 2:n-1
        dTdt(i) = alpha*(-(T(i) - T(i-1))/ dx^2 + (T(i+1)- T(i))/ dx^2);
    end
    dTdt(1) = alpha*(-(T(1) - T1s)/ dx^2 + (T(2)- T(1))/ dx^2);
    dTdt(n) = alpha*(-(T(n) - T(n-1))/ dx^2 + (T2s- T(n))/ dx^2);
    T = T+dTdt*dt;

    figure(1)
    plot(x,T,'Linewidth', 3)
    axis([0 L 0 50])
    xlabel('Distance ')
    ylabel('Temperature(\circC)')
    pause(0.1)
end

```

## Crank-Nicolsan

```python
# matlab
% heat with crank

clear;

cond = 1/2;   % heat conductivity
L = 1;    % length of the wire
T = 1;  %final time

Nt = 2500; % number of time steps
Dt = T/Nt;   % Time step
Nx = 50;   % number of space steps
Dx = L/Nx;    % space step
b = cond*Dt/(Dx*Dx);     % alpha parameter in the fd impleme

% the initial condition the initial temperature is of the sine function

for n = 1:Nx+1
    x(n) = (n-1)*Dx;
    u(n,1) = sin(pi*x(n));

end

% boundary conditions
for j = 1:Nt+1   % vector t
    u(1,j) = 0;
    u(Nx+1, j) = 0;
    t(j) = (j-1)*Dt;  % vector t due to time discrete

end

% Defining the M_right and M_left matrices in the crank nic methd
% as Mr and Ml are tridiagonal matrices , we first define the three
% diagnols and then build the matrices
aal(1:Nx-2) = -b;
bbl(1:Nx-1) = 2.+2.*b;
ccl(1:Nx-2) = -b;
MMl = diag(bbl,0) + diag(aal, -1)+ diag(ccl,1);

aar(1:Nx-2) = b;
bbr(1:Nx-1) = 2.-2.*b;
ccr(1:Nx-2) = b;
MMr = diag(bbr,0) + diag(aar, -1)+ diag(ccr,1);

% Implementing the method
for j = 1:Nt
    u(2:Nx, j+1) = inv(MMl)*MMr*u(2:Nx,j);
end

% graph
figure(1)
plot(x,u(:,1), '-b', x,u(:,round(Nt/100)),'--g',x,u(:,round(Nt/10)), ':b',x,u(:,Nt),'-.r')

title("temperature using crank")

xlabel('X')
ylabel('T')

% 3d
figure(2)
mesh(t,x,u)

title('Temperature with crank')
xlabel('t')
ylabel('x')

```

## Laplace

```python
%%%% define the rectangle and grid parameters
Lx = 1; % rectangle dimension
Ly = 1;
Ny = 100;
Nx = 100; % of gridpoints in x,y direction
ny = Ny+1;
nx = Nx+1;
dx = Lx/Nx;
dy = Ly/Ny;
x = (0:Nx)*dx;
y = (0:Ny)*dy;

boundary_index= [      1:nx,       1:nx:1+(ny-1)*nx, ....
    1+(ny-1)*nx:nx*ny,     nx:nx:nx*ny                ];

diagonals = [4*ones(nx*ny,1), -ones(nx*ny,4)];
A = spdiags(diagonals, [0 -1  1 -nx  nx], nx*ny, nx*ny);
I = speye(nx*ny);
A(boundary_index,:) = I(boundary_index,:);

%% set up boundary
b= zeros(nx, ny);
b(:,1) = 0;
b(1,:) = 0;
b(:,ny) = 4*x.*(1-x);
b(nx,:) = 0;
b = reshape(b, nx*ny,1);

Phi = A\b;
Phi = reshape(Phi, nx,ny);

[X,Y] = meshgrid(x,y);

v = [0.8  0.6  0.4  0.2  0.1  0.05  0.01];
contour(X,Y, Phi',v,'ShowText','on');
axis equal;

set(gca, 'YTick', [0  0.2  0.4  0.6  0.8  1]);
set(gca,  'XTick', [0  0.2  0.4   0.6  0.8  1]);

xlabel('$x$', 'Interpreter', 'latex','FontSize', 14);
ylabel('$y$', 'Interpreter', 'latex','FontSize', 14);
title('Solution of the Laplace equation','Interpreter', 'latex', 'FontSize',16);

```


text:
# Circuit Calculations
---

> ### Imports

```
import math
import sympy
from sympy.abc import *
pretty = sympy.pretty
sqrt = math.sqrt
```

> ### Resonant Frequency
* ω = 1/sqrt(L*C)
    * L = inductance, C = capacitance

```
L = 1 #@param {type:"number"}
C = 1 #@param {type:"number"}
w = 1/sqrt(L*C)
print(pretty(omega) + ' = ' + str(w))
```

    ω = 1.0

ext: ipynb


text:
Our solver will solve the Navier Stokes equations in two dimensions to calculate the velocities at each grid point and then move to the next time step. The continuity equation will be used as a governing equation for pressure by some manipulation. Hence the we will solve the following equations in theory:
$$
\begin{equation}
    \cfrac{\partial u}{\partial x} +\cfrac{\partial v}{\partial y}=0
\end{equation}
\\
\begin{equation}
    \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)
\end{equation}\\
\begin{equation}
    \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)
\end{equation}
$$
But since we do not have an explicit governing equation for pressure, we will take the divergence of the momentum equations, add them and use the continuity equation to eliminate the $\cfrac{\partial u}{\partial x} +\cfrac{\partial v}{\partial y}$ terms in the equations, to get the pressure poisson equation, given by,
$$
\begin{equation}
    \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)
\end{equation}
$$
This equation essentially ensures that continuity is satisfied at the current time step and also couples the pressure and velocities together.

## Constants and initial conditions

```python
import numpy as np
from matplotlib import pyplot as plt, cm
from solver import *

lx = 5      #length in the x direction
ly = 5      #length in the y direction
nx = 51     #grid points in x direction
ny = 51     #grid points in y direction
t = 10       #final time to calculate to
#nit = 10    #pseudo-time steps in each timestep
dx = lx / (nx - 1)
dy = ly / (ny - 1)
x = np.linspace(0, lx, nx)
y = np.linspace(0, ly, ny)
X, Y = np.meshgrid(x, y)

rho = 1
nu = 0.5
dt = .001 #time step size

nt=int(t/dt)
ds=int(100)

u = np.zeros((ny, nx)) #initial conditions
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
#b = np.zeros((ny, nx))
```

## Function for calculations of pressure poisson equation

```python
u = np.zeros((ny, nx)) #comment out these lines if you want to run from last time step
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))

vel_mag=(u**2+v**2)**0.5

np.savetxt('p0.csv', p, delimiter=',')
np.savetxt('u0.csv', u, delimiter=',')
np.savetxt('v0.csv', v, delimiter=',')
np.savetxt('Umag0.csv', vel_mag, delimiter=',')

#PRESSURE
pLeft=['N',0]
pRight=['N',0]
pTop=['D',0]
pBottom=['N',0]

pBCs=[pLeft,pRight,pTop,pBottom]
#pressureBC(pLeft,pRight,pTop,pBottom,p,dx,dy)

#VELOCITY

#x-component
uLeft=['D',0]
uRight=['D',0]
uTop=['D',10]
uBottom=['D',0]

uBCs=[uLeft,uRight,uTop,uBottom]
#XVelBC(uLeft,uRight,uTop,uBottom,u,dx,dy)

#y-component
vLeft=['D',0]
vRight=['D',0]
vTop=['D',0]
vBottom=['D',0]

#YVelBC(vLeft,vRight,vTop,vBottom,v,dx,dy)
vBCs=[vLeft,vRight,vTop,vBottom]

# solving
u, v, p = flow_solver(nt, u, v, dt, ds, dx, dy, p, rho, nu, pBCs,uBCs,vBCs)

```

```python
clearResults()
```

## Plotting the results

```python
fig = plt.figure(figsize=(11,7),dpi=100)
# plotting the pressure field as a contour
plt.contourf(X, Y, u, alpha=0.8, cmap=cm.viridis)
plt.colorbar()
## plotting the pressure field outlines
#plt.contour(X, Y, p, cmap=cm.viridis)
## plotting velocity field
plt.quiver(X[::lx, ::ly], Y[::lx, ::ly], u[::lx, ::ly], v[::lx, ::ly])
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.gca()
ax.set_aspect(1)
plt.show()
```

```python
fig = plt.figure(figsize=(11,7),dpi=100)
# plotting the pressure field as a contour
#plt.contourf(X, Y, u, alpha=0.8, cmap=cm.viridis)
## plotting the pressure field outlines
plt.pcolormesh(X, Y, u, alpha=1, cmap=cm.viridis)
plt.colorbar()
## plotting velocity field
#plt.streamplot(X[::lx, ::2*ly], Y[::lx, ::2*ly], u[::lx, ::2*ly], v[::lx, ::2*ly])
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.gca()
ax.set_aspect(1)
plt.show()
```

```python
import pandas as pd
import matplotlib.animation as animation

fig, ax = plt.subplots(dpi=100)
#levels=np.linspace(0,1.5,10)
def animate(i):
    tmpName='Umag'+str(i*ds)+'.csv'
    Tfile=pd.read_csv(r"Results/"+str(tmpName),header=None)
    Tarr=Tfile.to_numpy()
    # plotting the temperature field
    re = plt.contourf(X, Y, Tarr, 20, alpha=0.8, cmap=cm.viridis)
    #plt.xlabel('X')
    #plt.ylabel('Y')
    ax = plt.gca()
    ax.set_aspect(1)
    return re

anim = animation.FuncAnimation(fig, animate, frames=nt//ds, repeat=False)
anim.save(r'Gifs/cavity_Umag.gif', writer='imagemagick', fps=10)

```