left=[0,20,40,60,80,0];
top=[0,50,20,10,0,0];
bot=[0,0,0,0,0,0];
right=left;
dens=2;ord=6;

botl=length(bot); rightl=length(right);
topl=length(top); leftl=length(left);
len=botl;        % number of points per original row/column
m=(len-1)*dens;  % number of subintervals per row/column
n=m-1;           % number of unknowns per row/column
l=m+1;           % number of points per row/column

if dens > 1
  cobot=polyfit(1:len,bot,ord);
  coright=polyfit(1:len,right,ord);
  cotop=polyfit(1:len,top,ord);
  coleft=polyfit(1:len,left,ord);

  newpoints=(((1:l) - 1)/dens) + 1;    % preserve endpoints as 1 & len

  bot=polyval(cobot,newpoints);
  right=polyval(coright,newpoints);
  top=polyval(cotop,newpoints);
  left=polyval(coleft,newpoints);
end;

u=lapsolve1(bot,right,top,left);

surf(u);
axis ij
%title('temperature distribution')
