I have a system of coupled first order partial differential equations. Minimal code is shown below:

```
de1 = D(f(t, p), t) ==
p/t D( 1/3 f(t, p) + 2/15 g(t, p), p) + 2/5 1/t g(t, p);
de2 = D(g(t, p), t) ==
p/t D(11/21 g(t, p) + 2/3 f(t, p), p) + 2/7 1/t g(t, p) - g(t, p);
IC = {f(t0, p) == Exp(-p), g(t0, p) == 1/2 Exp(-p)};
t0 = 0.1;
sol = NDSolve({de1, de2, IC}, {f, g}, {t, t0, 1}, {p, 10^-2, 10},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}) // Flatten
```

First, note that the initial conditions are only specified on one boundary.

I want to export the values of functions f and g at specific t (say t1) in format {p, f(t1,p), g(t1,p)}. Also, is it possible to export at any t1 (not just at the evaluated grid)?

I am newbie in solving PDE’s. Any suggestions to improve the code is much appreciated.