initial_profile.m

function u=initial_profile(X,Y,M,N)
  IC=3;
  switch IC
    case 1,
      Z=0.025.*exp(-25.*(X-0.3).^2-25.*(Y-0.3).^2)+...
        0.03.*exp(-50.*(X-0.55).^2-50.*(Y-0.25).^2)+...
        0.035.*exp(-125.*(X-0.75).^2-125.*(Y-0.65).^2)+...
        0.02.*exp(-100.*(X-0.45).^2-100.*(Y-0.8).^2)+...
        0.025.*exp(-50.*(X-0.55).^2-50.*(Y-0.55).^2);
    case 2,
      Z=0.5.*exp(-100.*(X-0.25).^2-100.*(Y-0.3).^2)+...
        0.65.*exp(-50.*(X-0.5).^2-50.*(Y-0.25).^2)+...
        0.55.*exp(-30.*(X-0.7).^2-20.*(Y-0.65).^2)+...
        0.6.*exp(-80.*(X-0.4).^2-80.*(Y-0.8).^2)+...
        0.4.*exp(-90.*(X-0.5).^2-90.*(Y-0.55).^2)+...
        0.45.*exp(-100.*(X-0.8).^2-100.*(Y-0.3).^2);
    case 3,
      Z=0.3.*exp(-200.*(X-0.2).^2-200.*(Y-0.2).^2)+...
        0.35.*exp(-300.*(X-0.8).^2-300.*(Y-0.8).^2)+...
        0.4.*exp(-600.*(X-0.3).^2-600.*(Y-0.7).^2)+...
        0.25.*exp(-400.*(X-0.7).^2-400.*(Y-0.2).^2);
    case 4,
      Z=zeros(size(X));
      for i=1:10
        Z(2,floor(N*rand+1))=0.25*rand;
      end
  end
  u=[ones(M*N,1);reshape(Z,M*N,1)];
return