# Nonlinear System with Cross-Coupling Between Components

This example shows how to solve a nonlinear PDE system of two equations with cross-coupling between the two components. The system is a Schnakenberg system

`$\begin{array}{l}\frac{\partial {\mathit{u}}_{1}}{\partial \mathit{t}}-{\mathit{D}}_{1}\Delta {\mathit{u}}_{1}=\kappa \left(\mathit{a}-{\mathit{u}}_{1}+{{\mathit{u}}_{1}}^{2}{\mathit{u}}_{2}\right)\\ \frac{\partial {\mathit{u}}_{2}}{\partial \mathit{t}}-{\mathit{D}}_{2}\Delta {\mathit{u}}_{2}=\kappa \left(\mathit{b}-{{\mathit{u}}_{1}}^{2}{\mathit{u}}_{2}\right)\end{array}$`

with the steady-state solution ${\mathit{u}}_{1\mathit{S}}=\mathit{a}+\mathit{b}$ and ${\mathit{u}}_{2\mathit{S}}=\frac{\mathit{b}}{{\left(\mathit{a}+\mathit{b}\right)}^{2}}$. The initial conditions are a small perturbation of the steady-state solution.

### Solution for First Time Span

First, create a PDE model for a system of two equations.

`model = createpde(2);`

Create a cubic geometry and assign it to the model.

```gm = multicuboid(1,1,1); model.Geometry = gm;```

Generate the mesh using the linear geometric order to save memory.

`generateMesh(model,"GeometricOrder","linear");`

Define the parameters of the system.

```D1 = 0.05; D2 = 1; kappa = 100; a = 0.2; b = 0.8;```

Based on these parameters, specify the PDE coefficients in the toolbox format.

```d = [1;1]; c = [D1;D2]; f = @(region,state) [kappa*(a - state.u(1,:) + ... state.u(1,:).^2.*state.u(2,:)); kappa*(b - state.u(1,:).^2.*state.u(2,:)) ]; specifyCoefficients(model,"m",0,"d",d,"c",c,"a",0,"f",f);```

Set the initial conditions. The first component is a small perturbation of the steady-state solution ${\mathit{u}}_{1\mathit{S}}=\mathit{a}+\mathit{b}$. The second component is the steady-state solution ${\mathit{u}}_{2\mathit{S}}=\frac{\mathit{b}}{{\left(\mathit{a}+\mathit{b}\right)}^{2}}$.

```icFcn = @(region) [a + b + 10^(-3)*exp(-100*((region.x - 1/3).^2 ... + (region.y - 1/2).^2)); ... (b/(a + b)^2)*ones(size(region.x))]; setInitialConditions(model,icFcn);```

Solve the system for times 0 through 2 seconds.

```tlist = linspace(0,2,10); results = solvepde(model,tlist);```

Plot the first component of the solution at the last time step.

`pdeplot3D(model,"ColorMapData",results.NodalSolution(:,1,end));`

### Initial Condition for Second Time Span Based on Previous Solution

Now, resume the analysis and solve the problem for times from 2 to 5 seconds. Reduce the magnitude of the previously obtained solution for time 2 seconds to 10% of the original value.

```u2 = results.NodalSolution(:,:,end); newResults = createPDEResults(model,u2(:)*0.1);```

Use `newResults` as the initial condition for further analysis.

`setInitialConditions(model,newResults);`

Solve the system for times 2 through 5 seconds.

```tlist = linspace(2,5,10); results25 = solvepde(model,tlist);```

Plot the first component of the solution at the last time step.

```figure pdeplot3D(model,"ColorMapData",results25.NodalSolution(:,1,end));```

Alternatively, you can write a function that uses the results returned by the solver and computes the initial conditions based on the results of the previous analysis.

`NewIC = @(location) computeNewIC(results)`
```NewIC = function_handle with value: @(location)computeNewIC(results) ```

Remove the previous initial conditions.

`delete(model.InitialConditions);`

Set the initial conditions using the function `NewIC`.

`setInitialConditions(model,NewIC)`
```ans = GeometricInitialConditions with properties: RegionType: 'cell' RegionID: 1 InitialValue: @(location)computeNewIC(results) InitialDerivative: [] ```

Solve the system for times 2 through 5 seconds.

`results25f = solvepde(model,tlist);`

Plot the first component of the solution at the last time step.

```figure pdeplot3D(model,"ColorMapData",results25f.NodalSolution(:,1,end));```

### Function Computing Initial Conditions

```function newU0 = computeNewIC(resultsObject) newU0 = 0.1*resultsObject.NodalSolution(:,:,end).'; end```