function whs = mgmesh_wave_plot_initialise(g_bs, g_cg, g_mg, mgmesh, ll, edgecol, theta)
    x1m = min(mgmesh.P(1, :));
    x1M = max(mgmesh.P(1, :));
    x2m = min(mgmesh.P(2, :));
    x2M = max(mgmesh.P(2, :));
    pad = .05*max(x1M- x1m, x2M- x2m);

    F = mgmesh.F{ll};
    np = mgmesh.np(ll);
    P = [mgmesh.P(:, 1:np); 0*g_bs.'];

    m = max(max(abs(g_bs)), 1);
    m = max(max(abs(g_cg)), m);
    m = max(max(abs(g_mg)), m);

    cm_wave = gray(21);
    cm_diff = hot();

    ax_bs = subplot(2, 6, [3,4]);
    ph_bs = patch('Faces', F.', 'Vertices', P.', ...
                  'FaceVertexCData', g_bs,'FaceColor', 'interp',...
                  'EdgeColor', edgecol);
    axis 'equal';
    xlim([x1m-pad, x1M+pad]);
    ylim([x2m-pad, x2M+pad]);
    colormap(ax_bs, cm_wave);
    clim([-m, m]);
    colorbar;
    title("backslash (bs)");

    ax_cg = subplot(2, 6, [1,2]);
    ph_cg = patch('Faces', F.', 'Vertices', P.', ...
                  'FaceVertexCData', g_cg,'FaceColor', 'interp',...
                  'EdgeColor', edgecol);
    axis 'equal';
    xlim([x1m-pad, x1M+pad]);
    ylim([x2m-pad, x2M+pad]);
    colormap(ax_cg, cm_wave);
    clim([-m, m]);
    colorbar;
    title("nested cg");

    ax_mg = subplot(2, 6, [5,6]);
    ph_mg = patch('Faces', F.', 'Vertices', P.', ...
                  'FaceVertexCData', g_mg,'FaceColor', 'interp',...
                  'EdgeColor', edgecol);
    axis 'equal';
    xlim([x1m-pad, x1M+pad]);
    ylim([x2m-pad, x2M+pad]);
    colormap(ax_mg, cm_wave);
    clim([-m, m]);
    colorbar;
    title("nested vcycle");

    ax_cg_d = subplot(2, 6, [8,9]);
    ph_cg_d = patch('Faces', F.', 'Vertices', P.', ...
                  'FaceVertexCData', abs(g_cg-g_bs),'FaceColor', 'interp',...
                  'EdgeColor', edgecol);
    axis 'equal';
    xlim([x1m-pad, x1M+pad]);
    ylim([x2m-pad, x2M+pad]);
    colormap(ax_cg_d,cm_diff);
    colorbar;
    title("diff cg-bs");

    ax_mg_d = subplot(2, 6, [10,11]);
    ph_mg_d = patch('Faces', F.', 'Vertices', P.', ...
                  'FaceVertexCData', abs(g_mg-g_bs),'FaceColor', 'interp',...
                  'EdgeColor', edgecol);
    axis 'equal';
    xlim([x1m-pad, x1M+pad]);
    ylim([x2m-pad, x2M+pad]);
    colormap(ax_mg_d, cm_diff);
    colorbar;
    title("diff vcycle-bs");

    sgtitle(sprintf("solver test for \\theta-leapfrog, \\theta = %.4f", theta));

    whs.ph_bs = ph_bs;
    whs.ph_cg = ph_cg;
    whs.ph_mg = ph_mg;
    whs.ph_cg_d = ph_cg_d;
    whs.ph_mg_d = ph_mg_d;
    whs.ax_bs = ax_bs;
    whs.ax_cg = ax_cg;
    whs.ax_mg = ax_mg;
    whs.ax_cg_d = ax_cg_d;
    whs.ax_mg_d = ax_mg_d;
end