How do I change the poisson equation to solve for an alternative answer in Matlab-CodePudding

I'm trying to redo the Poisson equation so that ∂²Φ /∂x² = 2x² - 0.5x exp(x). Every time I try and input the right-hand side of the equation it comes up with a syntax error, any help would very much be appreciated.

% Solving of 1D Poisson equation
% using finite differences

% Clearing memory and figures
clear all
clf

% Define 1D numerical model
xsize=100; % Model size in horizontal direction, m
Nx=101; % Number of grid points
dx=xsize/(Nx-1); % Gridstep, m
x=0:dx:xsize; % Coordinates of grid points, m

% Defining global matrixes
L=zeros(Nx,Nx); % Koefficients in the left part
R=zeros(Nx,1); % Right parts of equations


% Composing global matrixes
% by going through all grid points
for j=1:1:Nx
    % Discriminating between BC-points and internal points
    if(j==1 || j==Nx)
        % BC-points
        L(j,j)=1; % Left part
        R(j)=0; % Right part
    else
        % Composing Poisson eq. d2T/dx2=1
        %
        % ---T(j-1)------T(j)------T(j 1)---- STENCIL
        %
        % 1/dx2*T(j-1) (-2/dx2)*T(j) 1/dx2*T(j 1)=1
        %
        % Left part of j-equation
        L(j,j-1)=1/dx^2; % T(j-1)
        L(j,j)=-2/dx^2; % T(j)
        L(j,j 1)=1/dx^2; % T(j 1)
        % Right part
        R(j)=2x;
    end
end

% Solve the matrix
S=L\R;

% Reload solutions
T=zeros(1,Nx); % Create array for temperature, K
for j=1:1:Nx
    T(j)=S(j);
end

% Visualize results
figure(1); clf
plot(x,T,'o r') % Plotting results

Again, I appreciate any help in this matter

CodePudding user response：

The problem is clear: "input the right-hand side of the equation it comes up with a syntax error", to correct the error you should properly define the right-hand side of the equation:

R(j)=2*x(j)^2-0.5*x(j)   exp(x(j));

Since it depends on x, you should point to each element j when calculating the coefficients in the for loop on the left and right hand sides of the equation.

Analytical solution of the PDE is:

Changing this line of code and plotting the discretized and analytical solutions produces

% Visualize results
figure(1); clf
plot(x,T,'o r');
hold on
plot(x , 0.166667*x.^4-0.0833333*x.^3-268811714181613573321934397213431784538112*x exp(x)-1)

which seems correct. Differences are due to discretization itself.