Here is the code for MP 3

%%%heat equation

clear all

close all

x = 80; % size of the rod/beam

M = 500;

X = linspace(0,x,M);

t = 5; % time

c = 0.35; % heat equation constant

dx = X(2)-X(1);

dt = 0.01;

T = 0:dt:t;

T = logspace(0,t,length(T));

const = 0.35;

M = length(X) % number of mesh points

N = length(T) % number of time steps

r = c*dt/(dx**2)

% initialize U

A = zeros(1,M);

I = zeros(1,M);

A(M/2) = 10;

I(M/2) = 10;

for i = -10:10

A(M/2-i) = 100;

I(M/2-i) = 100;

endfor

% set matrix

d = zeros(1,M)’+(1-2*r);

upper = zeros(1,M-1)’+r;

Q = diag(d) + diag(upper,1) + diag(upper, -1);

U = zeros(N,M);

for i = 1:N

f = A.^2./(I + 3) – const*A;

A = A*Q – f;

g = A.^2 – I;

I = I*Q*const + g;

U(i,:) = A;

plot(X,A)

pause(0.00001)

endfor

pause

surf(X,T,U);

pause

imagesc(X,T,U);