Why numerical solutions is not consistent with the exact solutions?

Why numerical solutions is not consistent with the exact solutions? Also the error is too large. How to modify the code to get consistent results? See the attached image for the equation of motions and their exact solutions.
clear; clc;
tic;
c1 = 0; % real constant
c2 = 0; % real constant
a = 1; % real constant
b=0.1; % two values b=0.1, b=0.01
lambda = -1./(2*(1 + 1i)); % complex or imaginary
N = 1000;
X = linspace(-10, 10, N);
T = linspace(0, 10, N);

% Calculate initial conditions
[x, ~] = meshgrid(X, 0); % t = 0 for initial conditions
%psi_init = A * exp(-1i * lambda * x);
%phi_init = B * exp(1i * lambda * x);
nq_init = -a*(b.^5*(1+(0 + c2).^2).^2 – 8*a*(0+c2).*(x+c1) – 4*a*b.^4*(0+c2).*(1+(0+c2).^2).*(x+c1) – 4*a*b.^2*(0+c2).*(x+c1).*(3+(0+c2).^2 + a.^2*(x+c1).^2) + b*(-3+(0+c2).^4 + 6*a.^2*(x+c1).^2 + a.^4*(x+c1).^4 + 2*(0+c2).^2.*(3 + a.^2*(x+c1).^2)) + 2*b.^3*(-1 + (0+c2).^4 + a.^2*(x+c1).^2 + (0+c2).^2*(4 + 3*a.^2*(x+c1).^2)));
dq_init = 1 + (0 + c2).^2 + b.^2 + (b*(0 + c2) + a*(x + c1)).^2;
q_init = (nq_init ./ (dq_init .^ 2));
ns_init = -(-3 + (0 + c2).^2 + b.^2*(-3 + (0 + c2).^2) + 4*1i*a*(x+c1) – 2*a*b*(0 + c2).*(x+c1) + a.^2*(x+c1).^2).*(exp(1i*(a*x + b*0)));
s_init = ns_init ./ dq_init;

initial_conditions = [q_init(:); s_init(:)];

t_span = T;

% Solve the differential equation numerically
options = odeset(‘RelTol’, 1e-4, ‘AbsTol’, 1e-8);
[t, u] = ode45(@(t, u) dydt(t, u, N), t_span, initial_conditions, options);

% Extract q and s from the solution
q = reshape(u(:, 1:N), [length(t), N]);
s = reshape(u(:, N+1:end), [length(t), N]);

% Calculate the exact solution over the mesh grid
[x, t] = meshgrid(X, T);
nq = -a*(b.^5*(1+(t + c2).^2).^2 – 8*a*(t+c2).*(x+c1) – 4*a*b.^4*(t+c2).*(1+(t+c2).^2).*(x+c1) – 4*a*b.^2*(t+c2).*(x+c1).*(3+(t+c2).^2 + a.^2*(x+c1).^2) + b*(-3+(t+c2).^4 + 6*a.^2*(x+c1).^2 + a.^4*(x+c1).^4 + 2*(t+c2).^2.*(3 + a.^2*(x+c1).^2)) + 2*b.^3*(-1 + (t+c2).^4 + a.^2*(x+c1).^2 + (t+c2).^2*(4 + 3*a.^2*(x+c1).^2)));
dq = 1 + (t + c2).^2 + b.^2 + (b*(t + c2) + a*(x + c1)).^2;
ns = -(-3 + (t + c2).^2 + b.^2*(-3 + (t + c2).^2) + 4*1i*a*(x+c1) – 2*a*b*(t + c2).*(x+c1) + a.^2*(x+c1).^2).*(exp(1i*(a*x + b*0)));
exact_q = (nq ./ (dq .^ 2));
exact_s = ns ./ dq;

% Compute the error
error_q = abs(q’ – exact_q);
error_s = abs(s’ – exact_s);

% Plotting the numerical solution
figure;
subplot(2, 1, 1);
surf(t, x, abs(q’));
title(‘Numerical Solution for |q_{n+1} – q_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘|q_{n+1} – q_n|’);
shading interp;

subplot(2, 1, 2);
surf(t, x, abs(s’));
title(‘Numerical Solution for |s_{n+1} – s_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘|s_{n+1} – s_n|’);
shading interp;

% 3D surface plot for exact solution
figure;
subplot(2, 1, 1);
surf(t, x, abs(exact_q’));
title(‘Exact Solution for |q_{n+1} – q_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘|q_{n+1} – q_n|’);
shading interp;

subplot(2, 1, 2);
surf(t, x, abs(exact_s’));
title(‘Exact Solution for |s_{n+1} – s_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘|s_{n+1} – s_n|’);
shading interp;

% Plot the error
figure;
subplot(2, 1, 1);
surf(t, x, error_q’);
title(‘Error for |q_{n+1} – q_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘Error’);
shading interp;

subplot(2, 1, 2);
surf(t, x, error_s’);
title(‘Error for |s_{n+1} – s_n|’);
xlabel(‘t’);
ylabel(‘x’);
zlabel(‘Error’);
shading interp;

% 2D error plot at specific times
figure;
times_to_plot = [0, 5, 10]; % Example times to plot error
for i = 1:length(times_to_plot)
time_idx = find(t_span >= times_to_plot(i), 1);
subplot(length(times_to_plot), 1, i);
plot(X, error_q(:, time_idx), ‘r’, X, error_s(:, time_idx), ‘b’);
title([‘Error at t = ‘, num2str(times_to_plot(i))]);
xlabel(‘x’);
ylabel(‘Error’);
legend(‘|q_{n+1} – q_n| Error’, ‘|s_{n+1} – s_n| Error’);
end

% Generate error table for fixed times from 0 to 10 in increments of 0.5
fixed_times = 0:0.5:10;
error_table_q = zeros(length(fixed_times), 1);
error_table_s = zeros(length(fixed_times), 1);

for i = 1:length(fixed_times)
time_idx = find(t_span >= fixed_times(i), 1);
error_table_q(i) = mean(error_q(:, time_idx));
error_table_s(i) = mean(error_s(:, time_idx));
end

% Display the error table
error_table = table(fixed_times’, error_table_q, error_table_s, …
‘VariableNames’, {‘Time’, ‘Error_q’, ‘Error_s’});
disp(error_table);

toc;

function dydt = dydt(~, u, N)
dydt = zeros(2*N, 1);
q = u(1:N);
s = u(N+1:2*N);

% Handle boundary conditions
for n = 3:N-2
dydt(n) = 1/2 * s(n) * conj(s(n)) – 1/2 * s(n+1) * conj(s(n+1));
dydt(N+n) = 1/2 * (q(n+1) – q(n)) * (s(n+1) – s(n));
end
endWhy numerical solutions is not consistent with the exact solutions? Also the error is too large. How to modify the code to get consistent results? See the attached image for the equation of motions and their exact solutions.
clear; clc;
tic;
c1 = 0; % real constant
c2 = 0; % real constant
a = 1; % real constant
b=0.1; % two values b=0.1, b=0.01
lambda = -1./(2*(1 + 1i)); % complex or imaginary
N = 1000;
X = linspace(-10, 10, N);
T = linspace(0, 10, N);