[R] overlay shaded area in base r plot

[ani jaya]

Dear R-Help,

I try to overlay two standard deviations from two means. How to change
the overlaid area between two shaded areas into a specific color or
maybe a soft color in a plot?

> dput(mean1[1:20])
c(122.194495954369, 118.955256282505, 115.540991140893, 113.116216840647,
111.24053267553, 109.827890459103, 108.523652505026, 107.033183023752,
105.259229707182, 103.124055471975, 100.877524901322, 98.6376579858244,
96.5156796979913, 94.6968535964952, 93.0502858942909, 91.5977750931902,
90.2651576455091, 89.0734094965085, 87.9562971122031, 86.9443758094494
)

> dput(mean2[1:20])
c(85.9733226064336, 85.1790888261485, 84.4436873347186, 83.8339820163413,
83.2527178882373, 82.7270110163027, 82.2807280513966, 81.9025314782845,
81.4641947900397, 81.0208948190769, 80.5647455355762, 80.0587583244909,
78.9354501839752, 78.4263183538841, 78.3688690670842, 77.9632850076303,
77.6888937972956, 77.5689173832678, 77.5860621564674, 77.6375257226031

> dput(sd1[1:20])
c(17.0473968859792, 19.7055838512895, 20.0280169279762, 19.5651973323109,
19.4292448372276, 19.0540881390057, 18.9295826565011, 18.7039744205801,
18.688319090752, 18.6946307822843, 18.7427879417687, 18.8069804986113,
18.8794754969972, 18.9830959093513, 19.0754455885966, 19.1360320444124,
19.0546404913351, 19.1280788210156, 19.2014747953522, 19.2456454651155
)

> dput(sd2[1:20])
c(19.2537538439847, 19.2462114145511, 19.2253213677253, 19.2476344684752,
19.2494964552016, 19.2871242318276, 19.2695982918065, 19.261879086708,
19.3410606442566, 19.3257171172684, 19.5029937941038, 19.5286520380856,
19.472587329424, 19.4815048434245, 19.3056947381611, 19.2193415347369,
19.1367354761276, 19.062832883228, 19.0123017669597, 18.9504334825052
)

what i've done:

plot(mean1[1:20],type="l", ylim=c(170,60),las=1,
     xlab="Data",ylab=expression(bold("Val")),cex.axis=1.2,font=2,
cex.lab=1.2,lwd=4)
polygon(c(1:20,20:1),c(mean1[1:20]+sd1[1:20],mean1[20:1]),col="lightblue")
polygon(c(1:20,20:1),c(mean1[1:20]-sd1[1:20],mean1[20:1]),col="lightblue")
polygon(c(1:20,20:1),c(mean2[1:20]+sd2[1:20],mean2[20:1]),col="lightyellow")
polygon(c(1:20,20:1),c(mean2[1:20]-sd2[1:20],mean2[20:1]),col="lightyellow")
lines(mean1,lty=1,lwd=2,col="blue")
lines(mean2,lty=1,lwd=2,col="yellow")

what i want is like in this source, but using base r plot if possible:

https://datascience.stackexchange.com/questions/30913/plots-with-shaded-standard-deviation

Any help and suggestion is very welcome!

Best,
Ani