One way using insets:
Plots gr code
using Plots, LaTeXStrings
uₜᵣᵤₑ(x,y) = x + 0.5y; uₚₗₒₜ(x,y) = 2x + y
x, y = 0:0.05:4, 0.05:0.1:0.95
surface(x, y, uₜᵣᵤₑ, c=cgrad([:red,:red]), camera=(60,30), cb=false,fa=0.8)
surface!(x, y, uₚₗₒₜ, c=cgrad([:blue,:blue]), zlabel=L"$u(t,x_0)$",fa=0.8)
title!(L"Exact and Approximated Surfaces $M=2$")
scatter!([0 0 0], [-1 NaN -1], lims=(0,1),
inset=(1,bbox(0.05,0.1,0.15,0.15)), subplot=2, msw=0, marker=:square,
legendfontsize=8, framestyle=:none, fg_color_legend=nothing, legend=:left,
color=[:red :white :blue], label=" "^4 .* [L"x + 0.5y" "" L"2x + y"]
)
