### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

I'm working on a transient 2D heat equation model and am having a few
problems with the boundary conditions for my 2D plate.
here's my code:

%% IMPLICIT CRANK NICOLSON METHOD FOR 2D HEAT EQUATION%%

clc;
clear all;

% define the constants for the problem
M = 25; % number of time steps
L = 1; % length and width of plate
k = 0.02; % time step
N = 5; % number of elements in plate (NxN)
h = L/(N-1); % grid spacing
x = [0:h:L]; % vector of x values
y = [0:h:L]; % vector of y values
%
[xx yy] = meshgrid(x,y);
%
% set constants a and c
a = (1 + (2*k)/(h^2));
c = (-k/(2*(h^2)));

% set the matrix for H
A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),
-1);
B = c*eye(N);
C = zeros(N);
%
% Create the matrix H, composed of sub-matrices A,B and C
Ap = eye(N);
Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);
Cp = ones(N) - Ap - Bp;
%
% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp
with copies of A, B, and C %
H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);
%
% % set up initial conditions for plate
a = 4; % the initial cond value for the plate
u = a*ones(N,N);
%
% set the boundary conditions for the plate
u(1,:) = 0; % b.c. for left hand end of plate
u(:,1) = 0; % b.c. for right hand end of plate
u(:,N) = 0; % b.c. for top hand end of plate
u(N,:) = 0; % b.c. for bottom hand end of plate
%
u = u(:);
for m = 1:M
u = H\u
axis manual; axis([0 L 0 L 0 a]);
drawnow
uu = reshape(u,N,N)
surf(xx,yy,uu), %colorbar
end

I think I need to put the b.c.'s inside the main loop for them to be
'remembered' throughout the iteration process. At the moment it
forgets them as soon as the iteration begins.

Any ideas would be most appreciated,

Geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Sure. At each iteration, the line u = H\u will replace the contents of u
where you first put the boundary conditions. Although
I am not sure if it is intended or not...

markus

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

a
N-1),
and
to
at a constant value. So that the edges of the plate are a constant
predetermined temperature.
What do you recommend i do to achieve this?

geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Do you really want the temperature to be 4 everywhere at time 0 but 0 at the edges? Assuming you do...

To make points stay the same temperature at all times, all you have to do is take the rows of H corresponding to those points, and
put a 1 on the diagonal and a zero everywhere else.

One more thing. The sum of every row of H should be 1. But in your case, this is only true for the rows corresponding to the
internal nodes. However, since you are specifying constant T BCs at all the edges, my point above will fix that for you.

-Paul

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

aul Skoczylas wrote:
a
N-1),
and
to

Ideally, my final model will enable the user to pick a point (or
points) on the grid to be a source or sink of heat. You then set an
initial temperature for the rest of the grid (an initial condition,
say the ambient temperature of the plate) and then observe the heat
transfer.

At this stage, if I can have all the boundaries at a
certain temperature that remains constant and an initial value for
the internal nodes then I'll be more than happy.

With the H matrix, wouldn't having the sum of each row of H equal to
1 mean nothing would happen during the iteration? I'm well confused
now!!
I can see what your saying about setting the H matrix so that certain
points are constant throughout the iteration, but will the
neighbouring points recognise this?

a trick somewhere!

Geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Having the sum of each row equal to one does NOT mean nothing happens. It just means nothing wierd happens! Nothing will happen to
a specific node if the value on the diagonal is 1, and everything else is zero--but that's what you want for your constant
temperature BCs.

Exception: when you implement heat sources or sinks, I think that the rows corresponding to those nodes might not have a sum of
one.

NOTE: I don't think you're using Crank-Nicolson. That would need an extra step, which you are missing. (I've never actually done
a C-N solution of this problem--I've always used straight implicit--but I just checked a text.) This appears to be a straight
implicit method. The coefficients in the matrix don't seem quite right to me, either, but I haven't fully examined them, either. I
don't like the way that thermal diffusivity is not listed, but assumed to be (I think) 1.0. Not a very useful program if you can't
specify the actual diffusivity.

Anyway, try it! Assuming your coefficients are correct (check them again later), set the rows in your H matrix corresponding to the
edges/corners so that the diagonal is 1 and the other elements are 0. Leave the other rows alone (and those rows do already have a
sum of 1--with 1.64 on the diagonal, and four -0.16 values elsewhere). Note that 16 of your 25 rows correspond to edges--there's
only 9 real calculations going on. And those real calculations are unaffected by the trivial nature of the other rows in the
matrix.

-Paul

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Use e.g. u_bc instead of u for b.c.'s.

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

aul Skoczylas wrote:
time 0
you
and
in
the
T
(or
set
equal

I think I've got a good idea of how to implement the b.c.'s now. The
alternative labelling will help to keep my confusion to a minimum.
The method is supposed to be crank nicholson. I borrowed the approach
<http://www.cs.berkeley.edu/~demmel/cs267/lecture17/lecture17.html>
not sure where i've tripped up then??!!
I haven't started looking at the coefficients yet. I've included the
diffusivity and density in my constant k. In my code it is labelled
incorrectly.

Thanks again paul and markus for all your help

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

That link says (this is copied/pasted, and might not look right, depending on the font):

k k
(1 + 2*---) * U(i,j,m+1) - ----- * ( U(i-1,j,m+1) + U(i+1,j,m+1) +
h^2 2*h^2
U(i,j-1,m+1) + U(i,j+1,m+1) )

k k
= (1 - 2*---) * U(i,j,m) + ----- * ( U(i-1,j,m) + U(i+1,j,m) +
h^2 2*h^2
U(i,j-1,m) + U(i,j+1,m) )

= b(i,j,m)

This looks correct for C-N if h =dx=dy, and k=diffusivity*timestep.

The RHS of the equation is b, which is a function of U. In your code, the RHS is just U, as in your line: u = H\u (which is what
you would use for straight implicit, but for that, your coefficients wouldn't be right). You need to generate that b vector from u,
and then calculate the next time step as u=H\b. This needs to be done at every time step.

I take issue with the statement in that link that C-N is "the simplest reasonable implicit method". Straight implicit is simpler,
and at least as reasonable as the explicit method they looked at first. (More so, since like C-N, it's unconditionally stable.)
The advantage of C-N is that it's a bit more accurate, but certainly more complicated. C-N is 2nd order accurate in space and time.
The plain explicit and implicit methods are 2nd order in space, but only first order in time.

-Paul

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

aul Skoczylas wrote:

So if i inserted a line of code with the equation for determining b
from u, this would solve my problem? or is there more to implementing
crank nicholson than that?
Would you recommend persuing the cn method or just using straight
forward implicit?
Also if i am to use straight implicit, in what way are my
coefficients off?
Finally here is my latest code with b.c.s and i.c. sorted out. I'm
currently trying to allow for the introduction of a hot or cold spot:

%% IMPLICIT CRANK NICOLSON METHOD %%

clc;
clear all;

% define the constants for the problem
M = 25; % number of time steps
L = 1; % length and width of plate
k = 0.01; % constant
N = 30 % number of elements in plate (NxN)
h = L/(N-1); % grid spacing
x = [0:h:L]; % vector of x values
y = [0:h:L]; % vector of y values
%
[xx yy] = meshgrid(x,y);
%
% set constants a and c
a = (1 + (2*k)/(h^2));
c = (-k/(2*(h^2)));
%
% set the matrix for H
A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),
-1);
B = c*eye(N);
C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C
Ap = eye(N);
Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);
Cp = ones(N) - Ap - Bp;
%
% Define the H matrix using KRON to "replace" 1's in Ap, Bp, and Cp
with copies of A, B, and C
H = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);
%
% % set up initial conditions for plate
a = 4; % the initial condition value for the plate
u = a*ones(N,N);
%
% set the boundary conditions for the plate
top = 4;
u(1,:) = top; % boundary condition for top of plate
left = 4;
u(:,1) = left; % boundary condition for left of plate
bottom = 4;
u(N,:) = bottom;% boundary condition for bottom
right = ;
u(:,N) = right; % boundary condition for right of plate
%
% set hot or cold spot
source_sink = 2;
k = 5; l = 10;
x_posn = k:l;
p = 5; q = 10;
y_posn = p:q;
u(y_posn,x_posn) = source_sink;
%
u = u(:);
%
for m = 1:M
u = H\u;
% replace the b.c. values
u(1:N) = left; % left b.c.
u(((2*N):N:N*N)) = bottom; % bottom b.c.
u((N+1):N:(N*N)) = top; % top b.c.
u((((N-1)*N)+1):(N^2)) = right; % right b.c.
% replace source/sink value
u((N*(k-1)+p):(N*(k-1)+q)) = source_sink;
u((N*(k)+p):(N*(k)+q)) = source_sink;
u((N*(k+1)+p):(N*(k+1)+q)) = source_sink;
u((N*(k+2)+p):(N*(k+2)+q)) = source_sink;
u((N*(k+3)+p):(N*(k+3)+q)) = source_sink;
u((N*(l-1)+p):(N*(l-1)+q)) = source_sink;

end

% create a surf plot of the temperature (z-axis) against grid
position
axis manual; axis([0 L 0 L 0 a]),% pause;
drawnow
uu = reshape(u,N,N)
surf(xx,yy,uu), colorbar

What do you think? The part for replacing source/sink values is a
little crude. I couldn't figure out how to define the set temperature
for each column of u and therefore just did each in turn.

Thanks again,

Geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Geoff Halstead" < XXXX@XXXXX.COM > wrote in message news: XXXX@XXXXX.COM ...

I think that would do it.

It depends what your goal is. If you want to do it as simply as possible, I'd forget C-N, and go straight implicit. But if you're
going to be making a code which is going to be running a huge simulation and computer time and accuracy are both very important,
than C-N is better. (But if computer time is very important, there's a whole lot you need to do, because this is currently very
inefficient.)

Where you now have (again, assuming k=diffusivity*timestep and h=dx=dy) for C-N
a = (1 + (2*k)/(h^2));
c = (-k/(2*(h^2)));
You would need for fully implicit
a = (1 + (4*k)/(h^2));
c = (-k/(h^2));

You really should try to understand where these come from! They are from a simple discretization of the heat conduction equation.
You simply insert numerical derivatives into the equation to replace the derivatives that are there. The difference between
explicit, implicit and C-N is at what time you evaluate the spatial derivatives.

Definitely crude! I really don't like the idea of replacing the temperatures at every step like that. But I've already made my
recommendation for how to do this instead... And since your sink/source is really only another constant temperature BC, the same
would work for it. Why make a change at every time step if you only need to do it once? Very inefficient!

I also get very concerned when you reuse variables. In your original code, you reused a, and that caused me some confusion. In
this, you still have that, but now you're also reusing k, and that is very dangerous--especially if you decide to continue with C-N,
since you're going to need the original k value!

Note that the improper implementation of C-N is causing errors. It looks to me like this plate is at 4 everywhere at the start.
The edges are held at 4 throughout, while a block in the middle is held at 2 after time 0. Therefore there should never be any
temperature on the plate which is not between 2 and 4. But when I run your code, the temperature is less than 2 everywhere except
in the regions where it is specified as 2 or 4.

Aside from all that, here's a way to speed it up:

The line u=H\u can be replaced with:
[LL,UU,PP]=lu(H);
u=UU\(LL\(PP*u));

(Of course, if you fix the C-N implementation, the u on the right will have to change to b.)

The trick here is that since H never changes, the first line can be put outside the loop, and only the second part goes inside the
loop. On my computer, doing this dropped the execution time by a factor of 10!

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

aul Skoczylas wrote:
determining
straight
N-1),
and
a
repeated ones. The new line for u = H\u really helps!
I'm now trying to change my code from crank nicholson to a simple
implicit scheme by altering the H matrix.
I'm going to impliment you idea for b.c.s later on, when I get the
code working and doing what I want it to.

cheers,

Geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

aul Skoczylas wrote:
determining
straight
N-1),
and
a
I was just wondering which implicit scheme you would recommend for me
to use. I eventually plan to translate this problem into 3D, so I'm
not sure whether this will have a bearing on the decision.
How easy is it to include the b(m) for c-nicholson in the code i have
at the moment?

cheers, geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

gt;> Definitely crude! I really don't like the idea of replacing
the
C-N,
when
right
can
a

just a quick one. I've altered my code so it is now crank nicholson
implicit. I'm now having a few problems with my constant k. When its
very small the results are what i expect but if i set k to anything
above 0.00001 it starts doing strange things on the corners of the
boundaries. what do you think is going wrong? here's the latest:

% define the constants for the problem
M = 100; % number of time steps
L = 1; % length and width of plate
k = 5; % constant
N = 30; % number of elements in plate (NxN)
h = L/(N-1); % grid spacing
x = [0:h:L]; % vector of x values
y = [0:h:L]; % vector of y values
%
[xx yy] = meshgrid(x,y);
%
% set constants a and c
a = (1 + (2*k)/(h^2));
c = (-k/(2*(h^2)));
d = (1 - (2*k)/(h^2));
e = (k/(2*(h^2)));
%
% set the matrix for H
A = a*eye(N) + c*diag(ones(1, N-1), 1) + c*diag(ones(1, N-1),
-1);
B = c*eye(N);
D = d*eye(N) + e*diag(ones(1, N-1), 1) + e*diag(ones(1, N-1),
-1);
E = e*eye(N);
C = zeros(N);

% Create the matrix H, composed of sub-matrices A,B and C
Ap = eye(N);
Bp = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);
Dp = eye(N);
Ep = diag(ones(1, N-1), 1) + diag(ones(1, N-1), -1);
Cp = ones(N) - Ap - Bp;
%
% Define the H_1 matrix using KRON to "replace" 1's in Ap, Bp, and Cp
with copies of A, B, and C
H_1 = kron(Ap, A) + kron(Bp, B) + kron(Cp, C);
% Define the H_2 matrix using KRON to "replace" 1's in Ap, Bp, and Cp
with copies of A, B, and C
H_2 = kron(Dp, D) + kron(Ep, E) + kron(Cp, C);
%
% set up initial conditions for plate
i_c = 4; % the initial cond value for the plate
u = i_c*ones(N,N);
%
% set the boundary conditions for the plate
top = 4;
u(1,:) = top; % bc for top of plate
left = 4;
u(:,1) = left; % bc for left side of plate
bottom = 4;
u(N,:) = bottom; % bc for bottom of plate
right = 4;
u(:,N) = right; % bc for right side of plate
%
% set hot or cold spot
source_sink = 2;
r = 5; s = 10;
x_posn = r:s;

p = 5; q = 10;
y_posn = p:q;
u(y_posn,x_posn) = source_sink;
%
u = u(:);
%
[SS_1,UU_1,PP_1] = lu(H_1);
%
for m = 1:M

b = H_2*u;
u = UU_1\(SS_1\(PP_1*b));
% replace the b.c. values
u(1:N) = left; % left b.c.
u(((2*N):N:N*N)) = bottom; % bottom b.c.
u((N+1):N:(N*N)) = top; % top b.c.
u((((N-1)*N)+1):(N^2)) = right; % right b.c.
% replace source/sink value
u((N*(r-1)+p):(N*(r-1)+q)) = source_sink;
u((N*(r)+p):(N*(r)+q)) = source_sink;
u((N*(r+1)+p):(N*(r+1)+q)) = source_sink;
u((N*(r+2)+p):(N*(r+2)+q)) = source_sink;
u((N*(r+3)+p):(N*(r+3)+q)) = source_sink;
u((N*(s-1)+p):(N*(s-1)+q)) = source_sink;
end

geoff

### MATLAB >> boundary condition problem - 2D heat equation using crank nicholson

Without looking at the code, here's a thought:

Your k constant is really thermal diffusivity multiplied by the time step. I don't know what material you're trying to model, but
the diffusivity of steel is on the order of 0.00002 ms or so, depending on the grade. Diffusivity of wood is 0.0000001 ms or
so.

So looking at steel, if your k is 0.00001, then your time step is 0.5 seconds. That's a pretty big time step. In wood, it's much
larger still!

Both C-N and straight implicit are unconditionally stable at all time step sizes--they don't blow up if your time step is too
large--but they aren't unconditionally accurate. You will get better accuracy with small time steps. You are just seeing
inaccuracies related to large time steps. These inaccuracies should settle out with time--you should approach the correct steady
state solution eventually. (But like I said, I didn't look at your code--there certainly could be other errors in there!)

I'd recommend making a simplified problem to which you know what the transient solution should be, and making sure your model is
correct with it, before you trust it on a complex problem. (Unfortunately, transient problems of this nature usually require a sum
of an infinite series to calculate analytically; often each term in the infinite series requires a solution to a transcendental
equation.)

-Paul

```I am looking to adapt a simple 1 dimensional heat quation code. The
code I am currently using makes use of dirichlet boundary conditions,
that is to say the end temperatures of the rod are set at a specific
temperature. I would however like to use neumann b.c.'s which
effectively model heat flow in or out of the end of the rod. I have
looked for a method of applying a b.c. which is time dependent but
cannotfigure out how to apply it to my code. Here is my matlab code,
which uses the 2nd order implicit crank nicholson method:

%% IMPLICIT CRANK NICOLSON METHOD %%

clc;
clear all;

% define the constants for the problem
L = 1; % length of bar
k = 0.5; % thermal conductivity of bar
N = 3; % number of elements in bar
J = 200; % number of time steps
dx = L/(N+1);
x = [0:dx:L]; % vector of x values
s = 0.6*k; % s is the step size

% set the matrix for u and the initial conditions
Id = eye(N+2); % sets up a zeros matrix of the desired
dimensions
S = ones(N+2); %
S = tril(S,1);
S = triu(S,-1);
S = S - 3*Id;
S(1,1) = 0; S(1,2) = 0;
S(N+2,N+1) = 0; S(N+2,N+2) = 0 ;
A = Id +s*S;
A = inv(Id - s*S)*A; % Crank nihcolson part
C =A'
%
% % set up initial conditions for bar
u = ones(1,N+2);
u = u.*0; %input('initial condition = ');

% set the boundary conditions for the bar
u(1) = 1; %input('left boundary condition = '); % boundary
condition for left hand end of bar
u(N+2) = 0 %input('right boundary condition = '); % boundary
condition for right hand end of bar

for j = 1:J
u = u*C;
drawnow;
plot(x,u);
end

Any ideas as to how to create time dependent b.c.s would be most
appreciated,

Geoff
```

```How to specific boundary condition for heat equation that at one side the B.C. takes the form of dU/dt=0? thanks!
```

```The following script (from Trefethen, Spectral Methods in
MATLAB) solves the 2nd order wave equation in 2 dimensions
(u_tt=u_xx+u_yy) using spectral methods, Fourier for x and
Chebyshev for y direction.
On its rectangular domain, the equation is subject to
Neumann boundary conditions along the sides ( u_y(x,+-1,t)
=0 ) and periodic boundary conditions at the ends ( u
(+3,y,t)= u(-3,y,t) ).

The calculation in BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1
Ny+1],2:Ny) is a (pre)solution of a system of equations (I
suppose the one coming up in vv([1 Ny+1],:) = BC*vv(2:Ny,:)
near the end) and Dy's are Chebyshev spectral
differentiation matrices in y direction. This BC is then
used near the end of the code in
vv([1 Ny+1],:) = BC*vv(2:Ny,:) where the Boundary
conditions are actually implemented. The first and last
rows of the solution vv (corresponding to y=+-1) are being
updated and the Neumann conditions are being used. I don't
seem to be able to see what is EXACTLY happending there.
How is vv([1 Ny+1],:) = BC*vv(2:Ny,:) related to u_y(x,+-
1,t)=0?

Thanks,
Omid.

PS. This script requires the function cheb() which I'm also
posting, in case anyone wanted to actually try the code.

% p37.m - 2D "wave tank" with Neumann BCs for |y|=1

% x variable in [-A,A], Fourier:
A = 3; Nx = 50; dx = 2*A/Nx; x = -A+dx*(1:Nx)';
D2x = (pi/A)^2*toeplitz([-1/(3*(dx/A)^2)-1/6 ...
.5*(-1).^(2:Nx)./sin((pi*dx/A)*(1:Nx-1)/2).^2]);

% y variable in [-1,1], Chebyshev:
Ny = 15; [Dy,y] = cheb(Ny); D2y = Dy^2;
BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny);

% Grid and initial data:
[xx,yy] = meshgrid(x,y);
vv = exp(-8*((xx+1.5).^2+yy.^2));
dt = 5/(Nx+Ny^2);
vvold = exp(-8*((xx+dt+1.5).^2+yy.^2));

% Time-stepping by leap frog formula:
clf, plotgap = round(2/dt); dt = 2/plotgap;
for n = 0:2*plotgap
t = n*dt;
if rem(n+.5,plotgap)<1
subplot(3,1,n/plotgap+1), mesh(xx,yy,vv), view(-10,60)
axis([-A A -1 1 -0.15 1]), colormap(1e-6*[1 1 1]);
text(-2.5,1,.5,['t = ' num2str(t)],'fontsize',18)
set(gca,'ztick',[]), grid off, drawnow
end
vvnew = 2*vv - vvold + dt^2*(vv*D2x +D2y*vv);
vvold = vv; vv = vvnew;
vv([1 Ny+1],:) = BC*vv(2:Ny,:);
end

% CHEB  compute D = differentiation matrix, x = Chebyshev
grid

function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal
entries
D  = D - diag(sum(D'));                 % diagonal entries
```

```Hi guys

I need a bit from help, have to do a program using Crank Nicholson  the equation of the heat for a plate. The program serious seemed to this one, that it is of a bar.

Can someone orientate me on since doing it?

Thank you

function MWD=calorcrank(m,N,L,c,T)
h=L/m;k=T/N;
lan=c^2*k/h^2;
%red espacial
for i=1:m-1
x(i)=i*h;
end
%matriz A (coeficientes del sistema que tenemos que resolver)
%es una matriz ftr tridiagonal
A=zeros(m-1);
A(1,1)=2+2*lan; A(1,2)=-lan;
A(m-1,m-2)=-lan; A(m-1,m-1)=2+2*lan;
for i=2:m-2
A(i,i-1)=-lan; A(i,i)=2+2*lan; A(i,i+1)=-lan;
end
%matriz B (coeficientes del sistema que tenemos que resolver)
%es una matriz ftr tridiagonal
B=zeros(m-1);
B(1,1)=2-2*lan; B(1,2)=lan;
B(m-1,m-2)=lan; B(m-1,m-1)=2-2*lan;
for i=2:m-2
B(i,i-1)=lan; B(i,i)=2-2*lan; B(i,i+1)=lan;
end
%Calculo de W0, le llamamos W
%equivalentemente W=ff(x)
for i=1:m-1
W(i)=ff(x(i));
end
W=W';
cont=1;
MW(cont,:)=W';
%Esquema iterativo
for j=1:N
%calculo del vector variable en j, bb
bb(1)=k*F(x(1),(j-1/2)*k)+lan*a(j*k)+lan*a((j-1)*k);
bb(m-1)=k*F(x(m-1),(j-1/2)*k)+lan*b(j*k)+lan*b((j-1)*k);
for i=2:m-2
bb(i)=k*F(x(i),(j-1)*k);
end
cc=B*W+bb'; %terminos independientes
%solucion del sistema de ecuaciones lineales AW=cc
%en este caso lo resolvemos por el metodo de eliminacion
%gaussiana que es el que tiene MATLAB por defecto W=A\cc
W=A\cc;
cont=cont+1;
MW(cont,:)=W';
end
t=0:k:T;
vi=a(t);
vf=b(t);
MWD(:,1)=vi';
MWD(:,m+1)=vf';
MWD(:,2:m)=MW;

En otros ficheros Mfiles:
>> edit a.m

function y=a(t)
y=exp(-t).*(1-t);

>> edit b.m

function y=b(t)
y=0;

>> edit F.m

function y=F(x,t)
y=exp(-t).*(2-t).*(x-1);

>> edit ff.m

function y=ff(x)
y=1-x+sin(pi.*x);
```