function [v,z,x,y] = stream(n,m,iterations)
%
% This function should compute the stream function for 
% stagnation flow around a corner
%
a = 0.; b = 1.; c = 0.; d = 1.; 
%iterations=100;
dx = (b-a)/(n-1);
dy = (d-c)/(m-1);
%initialize the solution
v=zeros(n,m);
for i = 1:n,
x(i) = a + (i-1)*dx;
end
for j = 1:m,
y(j) = c + (j-1)*dy;
end
% set horizontal boundary conditions
for i = 1:n,
% set the bottom boundary (floor)
v(i,1) = 0;
% set the top boundary (inflow)
v(i,m) = y(m)*x(i);
end
%
% set verical boundary conditions
for j = 1:m,
% set the left boundary (wall)
v(1,j) = 0;
% set the right boundary (outflow)
v(n,j) = x(n)*y(j);
end
%
% set the internal values
for k = 1:iterations,
%v;
for i = 2:n-1,
for j = 2:n-1,
v(i,j) = ( v(i-1,j) + v(i+1,j) + v(i,j-1) + v(i,j+1) )/4.; 
end
end
end
diff=0;
for i = 1:n,
for j = 1:m,
z(i,j) = x(i)*y(j);
dum = abs(z(i,j)-v(i,j));
if (dum > diff), diff=dum;
end
end
end
fprintf('max difference is %12.8f\n',diff);
% plot contours
contour(x,y,v,50);