Appendix C:
Source code for bkweuler.cpp

}
return (oldx);
}
 
LD f(LD xn, LD xnm1, LD deltax, vector S, vector U, vector dudx)
{
LD result;
result = xnm1 - xn + (deltax * xn * dot(S, dudx)) /
(normal(subtract(S, multiply(xn,U))) - dot(S, U) + xn);
return (result);
}
 
//--------------------------------------
//| Main program |
//--------------------------------------
main()
{
int num_points = 0, rho_initial_x_i;
LD x[max_points], rho[max_points], deltax, initial_condition, oldrho;
vector S[max_points], F[max_points], P[max_points], U[max_points];
char flag, junk[10];
file_input (num_points, deltax, initial_condition, x, F, S, U);
rho_initial_x_i = (num_points - 1) / 2;
rho[rho_initial_x_i] = initial_condition;
P[rho_initial_x_i].a = 0;
P[rho_initial_x_i].b = initial_condition;
 
for (int i = rho_initial_x_i + 1; i < num_points; i++)
{
if (i == rho_initial_x_i + 1)
oldrho = initial_condition + 0.0006;
rho[i] = secant(deltax, S[i], U[i], dudx(U, i, deltax), oldrho, rho[i-1]);
oldrho = rho[i-1];
P[i] = multiply(rho[i], U[i]);
}
 
for (i = rho_initial_x_i - 1; i >= 0; i--)
{
rho[i] = secant(0.0-deltax, S[i], U[i], dudx(U, i, deltax), rho[i+2], rho[i+1]);
P[i] = multiply(rho[i], U[i]);
}
cout << "\n\n\n\n\n\nDo you want file output? [Y/N] ";
cin.get(flag);
cin.getline(junk,10);
if ((flag == 'Y') || (flag == 'y'))
file_output(P, num_points);
return (0); //end main
}