matlab code for 2nd order differential equation with boundary conditions

(d^2 u(x))/〖dx〗^2 =(γu(x))/(1+αu(x))
at x=0,u=1
at x=1,du/dx=0
α=0.001,γ=1,2.5,5,10,100 plot u verses x at differnet γ values


% Parameters
alpha = 0.001;
gamma = [1, 2.5, 5, 10, 100];

% Define the differential equation
du2_dx2 = @(x, u) (gamma * u) / (1 + alpha * u);

% Define the range of x values
x_span = [0, 1];  % from x=0 to x=1

% Initial conditions at x=0
x0 = 0;
u0 = 1;


% Boundary conditions at x=1 (du/dx = 0)
boundary_condition = 0;

% Define function to compute residuals for boundary condition
boundary_residuals = @(u_end) (u_end(1) - boundary_condition);

% Solve the ODE with boundary condition
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[x, u] = ode45(@(x, u) [u(2); du2_dx2(x, u(1))], x_span, [u0, 0], options);

% Adjust u(1) to satisfy the boundary condition du/dx(1) = 0
u_end = u(end, :);
u_end(2) = u_end(2) + boundary_residuals(u_end);
[t, u] = ode45(@(x, u) [u(2); du2_dx2(x, u(1))], x_span, u_end, options);

% Plot u versus x
plot(t, u(:, 1), 'b-', 'LineWidth', 0.2);
xlabel('x');
ylabel('u(x)');
title('Plot of u(x) versus x');
axis([0 1 0 1]);  % Set axis limits for x and u

Answer 1

[Torsten]

You have a boundary value problem, not an initial value problem. ode45 is not suited in this case. Use bvp4c or bvp5c instead.

% Parameters
alpha = 0.001;
gamma = [1, 2.5, 5, 10, 100];
hold on
for i = 1:numel(gamma)
    fcn = @(x,y)[y(2);gamma(i)*y(1)/(1+alpha*y(1))];
    bc = @(ya,yb)[ya(1)-1;yb(2)];
    guess = @(x)[1;0];
    xmesh = linspace(0,1,20);
    solinit = bvpinit(xmesh, guess);
    sol = bvp4c(fcn, bc, solinit);
    plot(sol.x,sol.y(1,:))
end
hold off
grid on