Numercial Probelm: matrix is close to singlar or badly scaled

Hi all,

I have a problem in solving a pde function

ur*(partial delta/partial x)^2-ux*(partial^2) delta/(partial ^2)x+2*sqrt(ur*P)=b*sin(phi-delta)

The value of phi=0.5*exp(x-Ng/2)^2.(Ng=64)

I use the central difference for partial delta/ partial x=(delta(i+1)-delta(i-1))/(2*h) and (partial^2) delta/(partial ^2)x=delta(i+1)-2*delta(i)+delta(i-1)/h^2.

But matlab says Matrix is close to singular or badly scaled and cannot give the answer. So I modify the equation by add (epsilon/t+b)*delta(i) at both sides, then I get the figure and if I change the value of epsilon to 50, the amplitude of the figure decreases, just like epsilon is a damping factor.

My questions are as follows. Why by adding something at both sides, the equation can be solved? Theorectically the change of epsilon won't affact the result since it appears both sides. Why it plays a role of a damping factor?

Thank you!!

my codes is attached.

myfunction.m

     global Ng ux ur b P m t h epsilon x N H
      Ng=64; 
      ux=6; 
      ur=1; 
      b=1*ux;     P=0.01; 
        t=0.1; % Time Steps;
        h=0.1; % Spacial Spacing;
        epsilon=0.5; 
        x=linspace(0,Ng-h,Ng/h)';
        N=length(x);      A=zeros(N);
      for counter=2:N-1; % Center for Delta
         A(counter,counter)=2*ux/h^2+epsilon/t+b;
         A(counter,counter-1)=-ux/h^2;
         A(counter,counter+1)=-ux/h^2;
      end;
      A(1,1)=2*ux/h^2+epsilon/t+b;
      A(1,2)=-ux/h^2;
      A(1,N)=-ux/h^2;

      A(N,N)=2*ux/h^2+epsilon/t+b;
      A(N,N-1)=-ux/h^2;
      A(N,1)=-ux/h^2;

      H=inv(A);
      phi0=0.5*exp(-0.1*(x-(Ng)/2).^2); % Disturbance in Phi
      delta_guess=zeros(N,1);
      delta0=algpde(phi0,delta_guess); % Corresponding Delta
      plot(x,delta0)alg.m
    function delta=algpde(phi,delta_guess);

    global Ng ux ur b P m t h epsilon x N H A

    error=[];
    delta_old=delta_guess;
    for i=1:25;
    for counter=2:N-1;
       y(counter,1)=(epsilon/t+b)*delta_old(counter)-ur/(4*h^2)*(delta_old(counter+1)-delta_old(counter-1))^2+...
                    b*sin(phi(counter)-delta_old(counter))-...
                    2*sqrt(ur*P)/(2*h)*(delta_old(counter+1)-delta_old(counter-1));
    end;
    y(1,1)=(epsilon/t+b)*delta_old(1)-ur/(4*h^2)*(delta_old(2)-delta_old(N))^2+...
           b*sin(phi(1)-delta_old(1))-...
           2*sqrt(ur*P)/(2*h)*(delta_old(2)-delta_old(N));
    y(N,1)=(epsilon/t+b)*delta_old(N)-ur/(4*h^2)*(delta_old(1)-delta_old(N-1))^2+...
           b*sin(phi(N)-delta_old(N))-...
           2*sqrt(ur*P)/(2*h)*(delta_old(1)-delta_old(N-1));

    delta_new=H*y;
    error=[error norm(delta_new-delta_old)];
    delta_old=delta_new;
    end;
    delta=delta_new;
    error