function [U,xv] = heat_pwlin(z,u0,tv,xv)
% u = heat_pwlin(z,u0,tv)
% u = heat_pwlin(z,u0,tv,xv)
% find the solution of heat equation u_t = u_{xx} where initial value u0 is
%   continuous, piecewise linear, constant for x<=z(1) and x>=z(end)
% inputs:
%     z: nodes for initial data (increasing)
%    u0: values of u0 at nodes
%    tv: t-values for which we want to find solution
%    xv: (optional) x-values for which we want to find solution
% outputs:
%     u: array with values of solution at points xv, tv 
%        each column corresponds to one tv value
% Examples: 
%   heat_pwlin([0 1 2 3],[1 -2 2 -1],0:.02:1);
%   heat_pwlin([0 1 2 3],[1 -2 2 -1],0.3,0.7) % find u for t=.3,x=.7
%   x=-1:.1:1; heat_pwlin(x,1-x.^2,0:.02:1);  % u0(x)=1-x^2 for x in [-1,1], 0 otherwise

% needs colorcurves.m
tv(tv==0)=1e-15;   % use small number instead of zero to avoid division by zero
if nargin==3       % if no xv provided, make reasonable choice for xv:
  d = z(end)-z(1); h = min(diff(z));             % h is smallest interval length
  xv = linspace(z(1)-d,z(end)+d,1+3*10*d/h);     % spacing should be <= h/10
end
[X,T] = ndgrid(xv,tv);
G = @(x) x.*(erf(x/2)+1)/2 + exp(-x.^2/4)/pi^.5; % antiderivative of g(x)
U = u0(1)*ones(size(X));                   % U = u0(1)
c = diff(u0)./diff(z);                     % slopes
p = diff([0;c(:);0]);                      % jumps of slopes
for j=1:length(z)
  U = U + p(j)*T.^.5.*G(T.^-.5.*(X-z(j))); % add (jump of slope at z_j)*V(x-z_j,t)
end
colorcurves(X,T,U)

%    g(x) = erf(x/2)/2 + 1/2
%    antiderivative is G(x) = x*(erf(x/2)+1)/2 + exp(-x^2/4)/pi^.5
%    U(x,t) = g(t^-.5*x) solution of heat eq for u0=Heaviside fct
%    S(x,t) fundamental solution, partial deriv. of U wrt x
%    V(x,t) = t^.5*G(t^-.5*x) is antiderivative wrt x of U
%    (we omit t argument of S,U,V below)
% solution formula for integral over [a,b] where u0 is linear, and u0'=c=constant:
%     int_[a,b] S(x-y)*u0(y) dy 
%   = int_[a,b] U(x-y)*c dy  - U(x-b)*u0(b) + U(x-a)*u0(a)
%   = (-V(x-b)+V(x-a))*c - U(x-b)*u0(b) + U(x-a)*u0(a)
% Because of continuity of u0(x) all terms with U(x-a) and U(x-b) will cancel
% when we add contributions from all intervals together.