Solving coupled ODEBVP using MATLAB

I have a system which consists of the following coupled second order ODEBVP

d2a1/dx2 = C1*a1*a2 - C2*(1-a2^2);

d2a2/dx2 = C3*a1*a2 - C4*(1-a2^2);

The boundary conditions are

At x = 0, a1 = 1, da2/dx = 0;

At x = 1, a1 = C5, da2/dx = 0;

a1, a2 - Variables, C1-C5 - Constants

I tried using bvp4c but I am not able to solve it. Can someone please suggest me how can I solve my ststem and get plots for a1 and a2 wrt x. I have also tried using finite difference for solving the system but was not able to solve.

Answer 1

[Bjorn Gustavsson]

There's an old "proverb" going something like "when in danger when in doubt run around in circles scream and shoot".

I vaguely remember seeing solutions to how to handle these kinds of BVPs here, but cannot find them right now, you have a higher level of motivation and might be "luckier" on that fron. Your easiest approach might be to use the shooting method and convert the problem into an initial-value-problem where you search for the initial values of a2 and da2dx at 0 that gives you the correct values at x == 1. I've toted up something like this:

function err = err_bvpp(a_free,idx_free,a_targets,idx_targets,ode,C,a0,t_span)
% error-function for your type of BVP
a0(idx_free) = a_free;                            % Insert the free shooting-parameters for the initial condition
[~,a_tmp] = ode45(@(t,a) ode(t,a,C),t_span,a0);   % Integrate the IVP
err = sum((a_targets-a_tmp(end,idx_targets)).^2); % Calculate the distance from the end-target

Which seems to work OK for my very much toy-model problem:

% Defining your type of ODE:
dadx = @(t,a,C) [a(3);a(4);C(1)*a(1)*a(2)-C(2)*(1-a(2)^2);C(3)*a(1)*a(2)-C(4)*(1-a(2)^2)];
% Coefficient matrix:
C = [1 1 1 1];
% "time"-span:
x_span = [0 1];
% Initial condition:
a0 = [1 0.1,0.1,0]; % Only the first and last should be used in the ODE-integration!
idx_free = [2,3];   % Index to the components of a0 we're "aiming" our shooting over/in 
a_target = [0.3,0];   % target values [a1, dx2dx]
idx_target = [1 4]; % index of the target components
a_free0 = [1 1];
a_Que = fminsearch(@(a_free) err_bvpp(a_free,...
                                      idx_free,...
                                      a_target,...
                                      idx_target,...
                                      dadx,C,a0,x_span),a_free0);
A0 = a0;
A0(idx_free) = a_Que;
[ta,a] = ode45(@(t,a)dadx(t,a,[1 1 1 1]),linspace(0,1,201),A0);
plot(ta,a)
legend('a1','a2','da1dx','da2dx')
a(end,:)
a(end,idx_target) - a_target

If you aren't too unlucky you should have an easy time modifying the coefficient-array and target to get a solution to your problem.

HTH