%% 18.303 Matlab script for Assignment 1
%% Sept 2005 by Matthew J. Hancock

% Usage: 
%>> ps1('probs1')
%>> ps1('probs2')
%>> ps1('probs3')
%>> ps1('probs4b')
%>> ps1('probs4c')
% choose the one corresponding to the problem number you want to see

function varargout=ps1(varargin)

close all;

[varargout{1:nargout}]=feval(varargin{:});


%------------------------------
function probs1

x=0:0.1:2;
fx=1/2*(1-x);
figure(1),plot(x,fx,'k',x-2,fx,'k',x-4,fx,'k',x+2,fx,'k',[-3,3],[0,0],'k--',...
    -2:2:2,zeros(size(-2:2:2)),'k.',-2:2:2,0.5*ones(size(-2:2:2)),'ko',-2:2:2,-0.5*ones(size(-2:2:2)),'ko');
xlabel('x','fontname','times','fontsize',14);
set(gca,'fontname','times','fontsize',14,'xlim',[-3,3],'ylim',[-0.5,0.5]);
set(gca,'plotboxaspectratio',[1,0.3,1]);
%saveas(1,'sine.eps');

fx=abs(1/2*(1-x));
figure(2),plot(x,fx,'k',x-2,fx,'k',x-4,fx,'k',x+2,fx,'k');
xlabel('x','fontname','times','fontsize',14);
set(gca,'fontname','times','fontsize',14,'xlim',[-3,3],'ylim',[0,0.5]);
set(gca,'plotboxaspectratio',[1,0.3,1]);
%saveas(2,'cosine.eps');


x=0:0.01:1; t0=0;
u0=uxt(x,0);
u02=uxt(x,0.2);
u05=uxt(x,0.5);
u1=uxt(x,1);

figure(3),plot(x,u0,'k--',x,u02,'k',x,u05,'k-.',x,u1,'k');
xlabel('x','fontname','times','fontsize',14);
ylabel('u(x,t_0)','fontname','times','fontsize',14);

%------------------------------
function u=uxt(x,t)

u= (3*sin(pi*x)*exp(-pi^2*t)-sin(3*pi*x)*exp(-9*pi^2*t))/2;

%----------------------------------------
function probs2

close all;

x=linspace(0,1,128);
t=10.^linspace(-4,-1,100);
N=100;
EP=0.1;
mm=1:N;
nn=2*mm-1;
Bnn=(-1).^(mm+1).*sin( nn*pi*EP/2 )./( nn*pi*EP/2 );

fnt=20;

%% 3D plot
[XX,TT]=meshgrid(x(1:2:128),t);
u_u0 = 2*heatsolTypeI(XX,TT,N,nn,Bnn);
figure(1),mesh(XX,TT,u_u0);
set(gca,'fontname','times','fontsize',fnt,'view',[150,50],'Xdir','reverse','xlim',[0,1],'ylim',[0,0.1],'zlim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
zlabel('u(x,t)/u_0','fontname','times','fontsize',fnt); 
title('3D plot of u(x,t)','fontname','times','fontsize',fnt); 
%saveas(1,'uxt3D.eps');

% level curves
figure(2);
[cs,hdl]=contour(XX,TT,u_u0,[4,2,1,0.5,0.25]); 
clabel(cs,hdl,'fontsize',15,'rotation',0,'labelspacing',1024)
set(gca,'xlim',[0,1],'ylim',[0,0.1],'fontname','times','fontsize',fnt);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
title('Level curves u(x,t)=const','fontname','times','fontsize',fnt); 
%saveas(2,'uxtlevel.eps');

% 2D u-x plot
u_u0_ta = 2*heatsolTypeI(x,0.001,N,nn,Bnn);
u_u0_tb = 2*heatsolTypeI(x,0.01,N,nn,Bnn);
u_u0_tc = 2*heatsolTypeI(x,0.1,N,nn,Bnn);
figure(3),plot(x,u_u0_ta,'b',x,u_u0_tb,'r',x,u_u0_tc,'m'); hold on;
plot((1/2-EP/2)*[1,1],[0,1/EP],'k--',(1/2+EP/2)*[1,1],[0,1/EP],'k--',[1/2-EP/2,1/2+EP/2],[1/EP,1/EP],'k--');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,1],'ylim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('u(x,t_0)/u_0','fontname','times','fontsize',fnt); 
title('u-x plot','fontname','times','fontsize',fnt); 
%saveas(3,'ux.eps');

% 2D u-t plot
u_u0_xa = 2*heatsolTypeI(1/2,t,N,nn,Bnn);
u_u0_xb = 2*heatsolTypeI(0.4,t,N,nn,Bnn);
u_u0_xc = 2*heatsolTypeI(0.1,t,N,nn,Bnn);
figure(4),plot(t,u_u0_xa,'b',t,u_u0_xb,'r',t,u_u0_xc,'m');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,0.1],'ylim',[0,10]);
xlabel('t','fontname','times','fontsize',fnt); 
ylabel('u(x_0,t)/u_0','fontname','times','fontsize',fnt); 
title('u-t plot','fontname','times','fontsize',fnt); 
%saveas(4,'ut.eps');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function probs3

close all;

x=linspace(0,1,128);
t=10.^linspace(-4,-1,100);
N=100;
EP=0.1;
mm=1:N;
nn=2*mm;
Bnn=2*(-1).^mm.*sin( mm*pi*EP )./( mm*pi*EP );

fnt=20;

%% 3D plot
[XX,TT]=meshgrid(x(1:2:128),t);
u_u0 = 1+heatsolTypeII(XX,TT,N,nn,Bnn);
figure(1),mesh(XX,TT,u_u0);
set(gca,'fontname','times','fontsize',fnt,'view',[150,50],'Xdir','reverse','xlim',[0,1],'ylim',[0,0.1],'zlim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
zlabel('u(x,t)/u_0','fontname','times','fontsize',fnt); 
title('3D plot of u(x,t)','fontname','times','fontsize',fnt); 
%saveas(1,'P3uxt3D.eps');

% level curves
figure(2);
[cs,hdl]=contour(XX,TT,u_u0,[4,2,1,0.5,0.25]); 
clabel(cs,hdl,'fontsize',15,'rotation',0,'labelspacing',1024)
set(gca,'xlim',[0,1],'ylim',[0,0.1],'fontname','times','fontsize',fnt);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
title('Level curves u(x,t)=const','fontname','times','fontsize',fnt); 
%saveas(2,'P3uxtlevel.eps');

% 2D u-x plot
u_u0_ta = 1+heatsolTypeII(x,0.001,N,nn,Bnn);
u_u0_tb = 1+heatsolTypeII(x,0.03,N,nn,Bnn);
u_u0_tc = 1+heatsolTypeII(x,0.1,N,nn,Bnn);
figure(3),plot(x,u_u0_ta,'b',x,u_u0_tb,'r',x,u_u0_tc,'m'); hold on;
plot((1/2-EP/2)*[1,1],[0,1/EP],'k--',(1/2+EP/2)*[1,1],[0,1/EP],'k--',[1/2-EP/2,1/2+EP/2],[1/EP,1/EP],'k--');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,1],'ylim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('u(x,t_0)/u_0','fontname','times','fontsize',fnt); 
title('u-x plot','fontname','times','fontsize',fnt); 
%saveas(3,'P3ux.eps');

% 2D u-t plot
u_u0_xa = 1+heatsolTypeII(1/2,t,N,nn,Bnn);
u_u0_xb = 1+heatsolTypeII(0.4,t,N,nn,Bnn);
u_u0_xc = 1+heatsolTypeII(0.1,t,N,nn,Bnn);
figure(4),plot(t,u_u0_xa,'b',t,u_u0_xb,'r',t,u_u0_xc,'m');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,0.1],'ylim',[0,10]);
xlabel('t','fontname','times','fontsize',fnt); 
ylabel('u(x_0,t)/u_0','fontname','times','fontsize',fnt); 
title('u-t plot','fontname','times','fontsize',fnt); 
%saveas(4,'P3ut.eps');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function probs4b

close all;

x=linspace(0,1,128);
t=10.^linspace(-4,0,100);
N=100;
nn=1:N;
Bnn=4*(-1).^(nn+1)./(pi*(2*nn-1));

fnt=20;

%% 3D plot
[XX,TT]=meshgrid(x(1:2:128),t);
u_u0 = heatsolTypeMX(XX,TT,N,nn,Bnn);
figure(1),mesh(XX,TT,u_u0);
set(gca,'fontname','times','fontsize',fnt,'view',[167,52],'Xdir','reverse','xlim',[0,1],'ylim',[0,1],'zlim',[0,1]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
zlabel('u(x,t)/u_0','fontname','times','fontsize',fnt); 
title('3D plot of u(x,t)','fontname','times','fontsize',fnt); 
%saveas(1,'P4buxt3D.eps');

% level curves
figure(2);
[cs,hdl]=contour(XX,TT,u_u0,[0.99, 0.9, 0.5,0.25,0.1]); 
clabel(cs,hdl,'fontsize',15,'rotation',0,'labelspacing',1024)
set(gca,'xlim',[0,1],'ylim',[0,1],'fontname','times','fontsize',fnt);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
title('Level curves u(x,t)=const','fontname','times','fontsize',fnt); 
%saveas(2,'P4buxtlevel.eps');

% 2D u-x plot
u_u0_ta = heatsolTypeMX(x,0.001,N,nn,Bnn);
u_u0_tb = heatsolTypeMX(x,0.01,N,nn,Bnn);
u_u0_tc = heatsolTypeMX(x,0.1,N,nn,Bnn);
u_u0_td = heatsolTypeMX(x,0.5,N,nn,Bnn);
u_u0_te = heatsolTypeMX(x,1,N,nn,Bnn);
figure(3),plot(x,u_u0_ta,'b',x,u_u0_tb,'r',x,u_u0_tc,'m',x,u_u0_td,'m',x,u_u0_te,'m'); hold on;
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,1],'ylim',[0,1]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('u(x,t_0)/u_0','fontname','times','fontsize',fnt); 
title('u-x plot','fontname','times','fontsize',fnt); 
%saveas(3,'P4bux.eps');

% 2D u-t plot
u_u0_xa = heatsolTypeMX(0.9,t,N,nn,Bnn);
u_u0_xb = heatsolTypeMX(0.5,t,N,nn,Bnn);
u_u0_xc = heatsolTypeMX(0.01,t,N,nn,Bnn);
figure(4),plot(t,u_u0_xa,'b',t,u_u0_xb,'r',t,u_u0_xc,'m');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,1],'ylim',[0,1]);
xlabel('t','fontname','times','fontsize',fnt); 
ylabel('u(x_0,t)/u_0','fontname','times','fontsize',fnt); 
title('u-t plot','fontname','times','fontsize',fnt); 
%saveas(4,'P4but.eps');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function probs4c

close all;

x=linspace(0,1,128);
t=10.^linspace(-4,log10(0.5),100);
N=100;
EP=0.1;
nn=1:N;
Bnn=8./(EP*pi*(2*nn-1)).*cos( (2*nn-1)*pi/4 ).*sin( (2*nn-1)*pi*EP/4);

fnt=20;

%% 3D plot
[XX,TT]=meshgrid(x(1:2:128),t);
u_u0 = heatsolTypeMX(XX,TT,N,nn,Bnn);
figure(1),mesh(XX,TT,u_u0);
set(gca,'fontname','times','fontsize',fnt,'view',[150,50],'Xdir','reverse','xlim',[0,1],'ylim',[0,0.5],'zlim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
zlabel('u(x,t)/u_0','fontname','times','fontsize',fnt); 
title('3D plot of u(x,t)','fontname','times','fontsize',fnt); 
%saveas(1,'P4cuxt3D.eps');

% level curves
figure(2);
[cs,hdl]=contour(XX,TT,u_u0,[4,2,1,0.5,0.25]); 
clabel(cs,hdl,'fontsize',15,'rotation',0,'labelspacing',1024)
set(gca,'xlim',[0,1],'ylim',[0,0.5],'fontname','times','fontsize',fnt);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('t','fontname','times','fontsize',fnt); 
title('Level curves u(x,t)=const','fontname','times','fontsize',fnt); 
%saveas(2,'P4cuxtlevel.eps');

% 2D u-x plot
u_u0_ta = heatsolTypeMX(x,0.001,N,nn,Bnn);
u_u0_tb = heatsolTypeMX(x,0.03,N,nn,Bnn);
u_u0_tc = heatsolTypeMX(x,0.1,N,nn,Bnn);
u_u0_td = heatsolTypeMX(x,0.5,N,nn,Bnn);
figure(3),plot(x,u_u0_ta,'b',x,u_u0_tb,'r',x,u_u0_tc,'m',x,u_u0_td,'m'); hold on;
plot((1/2-EP/2)*[1,1],[0,1/EP],'k--',(1/2+EP/2)*[1,1],[0,1/EP],'k--',[1/2-EP/2,1/2+EP/2],[1/EP,1/EP],'k--');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,1],'ylim',[0,12]);
xlabel('x','fontname','times','fontsize',fnt); 
ylabel('u(x,t_0)/u_0','fontname','times','fontsize',fnt); 
title('u-x plot','fontname','times','fontsize',fnt); 
%saveas(3,'P4cux.eps');

% 2D u-t plot
u_u0_xa = heatsolTypeMX(1/2,t,N,nn,Bnn);
u_u0_xb = heatsolTypeMX(0.4,t,N,nn,Bnn);
u_u0_xc = heatsolTypeMX(0.1,t,N,nn,Bnn);
figure(4),plot(t,u_u0_xa,'b',t,u_u0_xb,'r',t,u_u0_xc,'m');
set(gca,'fontname','times','fontsize',fnt,'xlim',[0,0.5],'ylim',[0,10]);
xlabel('t','fontname','times','fontsize',fnt); 
ylabel('u(x_0,t)/u_0','fontname','times','fontsize',fnt); 
title('u-t plot','fontname','times','fontsize',fnt); 
%saveas(4,'P4cut.eps');


%--------------------------------------------------
function ff=heatsolTypeI(x,t,N,nn,Bn)

ff=0;
for m=1:N
    ff=ff + Bn(m)*sin(nn(m)*pi*x).*exp(-pi^2*nn(m)^2*t);
end

%--------------------------------------------------
function ff=heatsolTypeII(x,t,N,nn,Bn)

ff=0;
for m=1:N
    ff=ff + Bn(m)*cos(nn(m)*pi*x).*exp(-pi^2*nn(m)^2*t);
end

%--------------------------------------------------
function ff=heatsolTypeMX(x,t,N,nn,Bn)

ff=0;
for m=1:N
    ff=ff + Bn(m)*cos((2*nn(m)-1)/2*pi*x).*exp(-pi^2*(2*nn(m)-1)^2*t/4);
end