Faster Way of Plotting a Multivariate Integral with No Closed Form and Maximized

I have a double integral that looks like this.

$$J=\int_0^2\int_{-2}^2 \sin(y\cdot B_z+z\cdot B_y+\phi_y+\phi_0)\,dy\,dz$$

$$\phi_y$$ I am defining as different functions of y, e.g. $$\phi_y=0$$, $$\phi_y=.5y$$, $$\phi_y= \left\{\begin{array}{ll} 0 & y\leq 0 \\ 1 & y\gt 0 \\ \end{array} \right.$$

$$B_y$$ and $$B_z$$ are the axes I am plotting over, while $$y$$ and $$z$$ are the real space coordinates of the surface I am integrating over.

So a simplified version of what I have written is

\*SubsuperscriptBox[$$\[Integral]$$, $$0$$, $$2$$]$$\*SubsuperscriptBox[\(\[Integral]$$, $$-2$$, $$2$$]Sin[\[Phi]0 +
y*Bz + z*By +  .2
\*SuperscriptBox[$$y$$, $$2$$]] \[DifferentialD]y \[DifferentialD]z\)\
\);
Plot[MaxValue[Limit[Ic, By -> 0], \[Phi]0], {Bz, -10, 10},
PlotRange -> All] // AbsoluteTiming


I let $$\phi_y$$ be 0 to have the computation be as simple as possible and only plot a slice of $$J$$ at a given $$B_y$$. And I am maximizing $$J$$ with respect to $$\phi_0$$. As you can see it takes a minute and a half to plot. This is completely fine, but if I try more complicated functions the computation can take much longer.

I need to try many different $$\phi_y$$ out, like $$\phi_y=\sin(y)$$, but these more complicated functions don't always have a closed form solution and might take an extremely long time to do. Additionally, I sometimes run into a problem where MaxValue gives errors that it does not evaluate to a numerical value, and I have not figured out exactly what causes those.

So my questions are..

1. What is a more efficient way of doing this process? I have thought about using NDSolve, approximations, or something similar, but haven't been able to get it to work yet.

2. Say I would like to tune $$\phi_y$$ with some parameter and watch the plot change using Manipulate or something similar. Is there a way to do that without running the computation over and over so I can see the evolution smoothly?

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 May 22 at 23:39
• People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this meta Q&A helpful – Michael E2 May 22 at 23:40
• My bad, I added a snippet of code to copy and paste. – user72882 May 23 at 2:53
• @user72882 Why do you use By->0? In that case we have 1D integral. Or it is just for simplification? Also don't use By since it is system symbol in Mathematica. Let put for instance B2 instead. – Alex Trounev May 23 at 23:35
• @AlexTrounev I am taking the limit as By->0 because it is easier for me to interpret what changing different parameters does. once I get a good grasp on that I will start visualizing the entire 2D surface. I did not realize I was typing the word 'By' either, thanks for pointing that out, I will change it to B2 in my code. – user72882 May 24 at 6:20