# Graphing a Custom Function

I'm trying to make a graph of a function involving several steps, steps that I can't put straight into the Plot[ ] function. Here is the outline of the process:

• Two real numbers, A and B, are defined, which determine the coefficients to a system of differential equations.

• Mathematica solves the system for me (initial conditions are supplied) and returns the solutions.

• I take one of the solutions and integrate it (a definite integral), which produces a number, the final output.

You can see this is a process which takes two inputs and produces an output, but not as a simple function. What I want to do is make a graph plotting the final output as a function of A and B. Apologies if there's a simple fix, I've tried to search the documentation thoroughly.

### Edit

Here is the code I am working with

A = 1.5; B = 2;

Comatrix =
{{-(A + 2*B), 3*B, 0, 0, 0, 0, 0, 0, 0},
{A, -(A + 3*B), 4*B, 0, 0, 0, 0, 0, 0},
{0, A, -(A + 4*B), 5*B, 0, 0, 0, 0, 0},
{0, 0, A, -(A + 5*B), 6*B, 0, 0, 0, 0},
{0, 0, 0, A, -(A + 6*B), 7*B, 0, 0, 0},
{0, 0, 0, 0, A, -(A + 7*B), 8*B, 0, 0},
{0, 0, 0, 0, 0, A, -(A + 8*B), 9*B, 0},
{0, 0, 0, 0, 0, 0, A, -(A + 9*B), 10*B},
{0, 0, 0, 0, 0, 0, 0, A, -(10*B)}};

Nmatrix[t_] = {N2[t], N3[t], N4[t], N5[t], N6[t], N7[t], N8[t], N9[t], N10[t]};

system = Derivative[1][Nmatrix][t] == Comatrix . Nmatrix[t];

sol =
DSolve[
{system,
N2[0] == 1, N3[0] == 0, N4[0] == 0, N5[0] == 0, N6[0] == 0, N7[0] == 0,
N8[0] == 0, N9[0] == 0, N10[0] == 0},
{N2, N3, N4, N5, N6, N7, N8, N9, N10}, t];

{N2ans[t_], N3ans[t_], N4ans[t_], N5ans[t_], N6ans[t_],
N7ans[t_], N8ans[t_], N9ans[t_], N10ans[t_]} =
{N2[t], N3[t], N4[t], N5[t], N6[t], N7[t], N8[t], N9[t], N10[t]} /. Flatten[sol];

f[t_] = 2*B*N2ans[t];
Integrate[t*f[t], {t, 0, Infinity}]

0.326222

As noted by J.M. in a comment, the solution to the ODEs is given by MatrixExp. So, the integrand in the final line of code is,

s = MatrixExp[Comatrix t] [[1,1]];

which is a rather lengthy RootSum. As verification,

NIntegrate[2 B s /. A -> 1.5 /. B -> 2], {t, 0, Infinity}] // AbsoluteTiming
(* {0.0733787, 0.326222} *)

which reproduces the sample result in about a tenth of a second. s can be used for plotting, for instance by

Flatten[ParallelTable[{A, B, NIntegrate[t 2 B s, {t, 0, Infinity}]},
{A, 0.1, 1.0, 0.1}, {B, 0.1, 1.0, 0.1}], 1];
ListPlot3D[%, PlotRange -> All, AxesLabel -> {A, B, f},
ImageSize -> Large, LabelStyle -> Directive[12, Bold, Black]]

• MatrixExp[] supports the "action form", so you can write the solution as First[MatrixExp[Comatrix t, UnitVector[9, 1]]]. Jul 26, 2017 at 18:01
• @J.M. Or, I could write MatrixExp[Comatrix t][[1,1]] Jul 26, 2017 at 18:07
• Yes indeed, but for a large enough matrix, the action form often takes less effort to compute since it only generates a vector and not the full exponential. Jul 26, 2017 at 18:09
• @J.M. Great point, which I did not know. Thanks. Jul 26, 2017 at 18:10
• One more question. I'm trying to graph A and B against the integral of 2Bt*N2, but the label on this graph is just N2. Is that a mislabelling? Jul 27, 2017 at 19:33