That is because the vectors for t and a are the same length in your particular example (they also have the same entries). If you change one of them, for example choose dt = tmax/1000.0 you will see the problem.
Yes, in heatmap you should first give the values for the x-axis and then the y-axis, but which is which has to match the actual data. In this case, sol_pde[S(t, a)] has the time points along the first dimension (which is “down”, so it corresponds to the y-axis when you would plot a matrix) and the space points along the second dimension (which is “sideways”, so it corresponds to the x-axis ).
So another solution to fix the plot would be to just transpose the matrix, i.e. plot transpose(soln) to flip the order of the axes without having to change anything else.
(I find it’s a little bit like USB-A plugs … I usually need at least two tries to figure out the right order of the axes when plotting heatmaps – choosing differently-sized ranges for the x- and y- directions can help to spot the difference )
Thanks for the hepl so far. With @mkretzschmar and @chvandorp I am now trying to extend the model to add a contact function c(a,b) to the model, specifying the rate at which individuals of age a make contact with individuals of age b. For this to work we would need to replace the integral over a in the above model (Ia(I(t,a))) with the integral (over the whole domain) over b of c(a,b)*I(b,t). So far we have been unsuccessful getting this to work.
)